Report an errorErrata
Part II · Module 11

Correlation & time-delay estimation

Speckle tracking, element-to-element aberration delays, sound-speed phase shifts and ULM bubble pairing all come down to one task: estimate a small shift between two noisy, partly decorrelated echoes. This module covers the estimators (normalized cross-correlation, subsample fits, phase of the complex correlation, Kasai), how biased they are, how small their variance can get (the Walker–Trahey bound), and where they fail (false peaks).

You will be able to

  • Explain why the cross-correlation peak is the ML delay estimate in white Gaussian noise, and when NCC is a usable quality metric
  • Choose a subsample interpolator and predict its bias curve (zero at 0 and 0.5 samples, maximal in between)
  • Evaluate the Walker–Trahey bound for given \(f_0\), \(B\), \(T\), SNR and \(\rho\), and recognise the false-peak threshold
  • Use lag-one autocorrelation (Kasai) on IQ, compute the Nyquist velocity, and reason about phase wrapping
  • Size a tracking kernel against strain decorrelation, and explain cumulative drift in neighbour-element phase-screen estimates

Assumes

  • Module 03 · Convolution, LTI systems and filters and Module 04: circular vs linear convolution, DFT, zero padding
  • Module 07 · Random vectors, covariance & correlation: correlation coefficient, covariance
  • Module 08 · Analytic signal & IQ demodulation: analytic signal, IQ demodulation, phase of a complex signal
  • Pulse-echo geometry: two-way delay \(\tau\) ↔ axial displacement \(\Delta z = c\tau/2\)
Exact formulaclosed form plotted as isSimulationMonte Carlo or estimator runs in this page (JS) or in the MATLAB/Octave scripts; solver named under the badgeApprox. modelsimplified physics (e.g. expected correlation of a Gaussian-spectrum speckle, a bound evaluated with an assumed decorrelation law); assumptions under the badgeSchematicdrawing for intuition, not computedSpeckle is synthesized in the frequency domain from a Gaussian pulse spectrum with seeded random scatterers. No Field II or measured channel data are used in the figures.

11.1Cross-correlation as a delay estimator

Model two echo segments from the same speckle region before and after a shift \(D\) (two-way time). With \(s(t)\) the common speckle signal and \(n_1\), \(n_2\) independent noise:

\[ x(t) = s(t) + n_1(t), \qquad y(t) = s(t - D) + n_2(t), \qquad t \in [0, T] \]
(11.1)
\(x(t),\ y(t)\)
reference and delayed records, V
\(s(t)\)
common speckle (or point-target) signal
\(n_1,\ n_2\)
independent electronic noise in the two records
\(D\)
true two-way delay of the second record, s
\(T\)
window (kernel) duration, s
\(t\)
time, s
The same model covers a displaced tissue window (elastography, ARFI), the signal on element \(i+1\) relative to element \(i\) (phase screen), and a bubble patch in frame \(k+1\) relative to frame \(k\) (ULM).

The cross-correlation function and its normalized form are

\[ R_{xy}(\tau) = \mathrm{E}\bigl[\, x(t)\, y(t + \tau) \,\bigr] = R_{ss}(\tau - D) \qquad \text{(noise independent ⇒ cross terms vanish)} \]
(11.2)
\(R_{xy}(\tau)\)
cross-correlation of \(x\) and \(y\) at lag \(\tau\)
\(\mathrm{E}[\cdot]\)
expectation over the noise and speckle ensemble
\(R_{ss}\)
autocorrelation of the common signal \(s\), maximal at zero lag
\(\tau\)
lag, s
\(x,\ y,\ D,\ t\)
records, delay and time, as in (11.1)
\[ \hat{\rho}(k) = \frac{\sum_{n \in W} x[n]\, y[n+k]}{\sqrt{\sum_{n \in W} x[n]^2 \;\cdot\; \sum_{n \in W} y[n+k]^2}}, \qquad \hat{D} = T_s \cdot \arg\max_k \hat{\rho}(k) \]
(11.3)
\(\hat{\rho}(k)\)
normalized cross-correlation (NCC) at lag \(k\), \(\lvert\hat\rho\rvert \le 1\)
\(x[n],\ y[n]\)
sampled reference and delayed records, as in (11.1)
\(W\)
set of sample indices in the kernel window
\(n\)
sample index
\(k\)
integer lag, samples
\(T_s\)
sampling interval \(1/f_s\), s
\(\hat{D}\)
delay estimate, s
\(\arg\max\)
lag at which \(\hat\rho\) is largest
The kernel \(x\) stays fixed; \(y\) slides. Because \(R_{ss}(\tau)\) is maximal at \(\tau = 0\) (Cauchy–Schwarz), the peak of \(R_{xy}\) sits at \(\tau = D\).

Why the peak is the maximum-likelihood estimate (sketch)

Treat \(s\) as known for a moment and only \(y\) as noisy (white Gaussian, variance \(\sigma^2\) per sample). The log-likelihood of the delay is

\[ \ln p(y; D) = -\frac{1}{2\sigma^2} \sum_n \bigl( y[n] - s[n - D] \bigr)^2 + \mathrm{const} \]
(11.4)
\(\ln p(y;D)\)
log-likelihood of the delay \(D\) given the record \(y\)
\(\sigma^2\)
variance of the white Gaussian noise per sample
\(y[n],\ s[n-D]\)
noisy record and delayed signal, as in (11.1)
\(n\)
sample index
\(\mathrm{const}\)
terms independent of \(D\)
\[ = -\frac{1}{2\sigma^2} \Bigl[ \sum y^2 + \sum s^2[n - D] - 2 \sum y[n]\, s[n - D] \Bigr] \]
(11.5)
\(\sigma^2,\ y,\ s,\ D\)
noise variance, record, signal and delay, as in (11.4)
\(\sum\)
sum over the sample index \(n\)
\(n\)
sample index
For a record longer than the signal, \(\sum s^2[n - D]\) does not depend on \(D\), so maximising the likelihood equals maximising the correlation \(\sum y[n]\, s[n - D]\): a matched filter evaluated over all candidate delays.

With \(s\) unknown and both records noisy, the ML estimator (Knapp & Carter, 1976) is a generalized cross-correlation: multiply the cross-spectrum \(G_{xy}(f)\) by the weight \(\lvert\gamma(f)\rvert^2 / \bigl[\lvert G_{xy}(f)\rvert\,(1 - \lvert\gamma(f)\rvert^2)\bigr]\) before the inverse transform, where \(\gamma\) is the coherence between \(x\) and \(y\). For white electronic noise and an in-band SNR well above 1, \(\lvert\gamma\rvert^2/(1 - \lvert\gamma\rvert^2) \approx S(f)/(2N_0)\) and \(\lvert G_{xy}\rvert = S(f)\), so the weight is nearly flat inside the transducer band and falls off outside it: plain correlation of band-pass-filtered RF is already close to ML. Filter out-of-band noise first; that is the pre-whitening step of the GCC.

NCC as a quality metric

For a long window, the expected NCC peak is the correlation coefficient between the two windows:

\[ \mathrm{E}\bigl[\hat{\rho}_{\max}\bigr] \approx \rho \cdot \frac{\mathrm{SNR}}{1 + \mathrm{SNR}} \]
(11.6)
\(\hat{\rho}_{\max}\)
peak of the NCC over lag
\(\rho\)
correlation of the noise-free signals (speckle decorrelation), 0 to 1
\(\mathrm{SNR}\)
signal-to-noise power ratio, equal in both records
\(\mathrm{E}[\cdot]\)
expectation
The second factor is the coherence loss due to noise. Thresholding the peak NCC (often 0.7–0.9) rejects unreliable estimates; Section 11.3 shows why that works: the jitter grows as \(\sqrt{1/\rho^2 - 1}\).

Biased vs unbiased correlation estimates, and FFT correlation

For a single record of length \(N\) the sample correlation can be normalized two ways:

\[ \begin{aligned} \hat{R}_b(k) &= \frac{1}{N} \sum_{n=0}^{N-1-\lvert k\rvert} x[n]\, y[n+k] \\ \hat{R}_u(k) &= \frac{1}{N - \lvert k\rvert} \sum_{n=0}^{N-1-\lvert k\rvert} x[n]\, y[n+k] \end{aligned} \]
(11.7)
\(\hat{R}_b(k),\ \hat{R}_u(k)\)
biased and unbiased sample cross-correlation at lag \(k\)
\(x[n],\ y[n]\)
records of length \(N\), as in (11.3)
\(N\)
record length, samples
\(n\)
sample index
\(k\)
lag, samples
\(\lvert k\rvert\)
absolute value of the lag
\(\mathrm{E}[\hat{R}_b(k)] = (1 - \lvert k\rvert/N)\, R(k)\): a triangular taper that pulls the peak toward \(k = 0\) (a delay bias that is small only when \(\lvert D\rvert \ll N\)). \(\hat{R}_u\) is unbiased but its variance blows up at large lags, and it is not guaranteed positive semidefinite. NCC with a fixed kernel and a sliding search window (as used in the demos) avoids both: every lag uses exactly \(W\) products.

FFT correlation computes \(X^*(f)\,Y(f)\) and inverts it. The DFT makes this a circular correlation (Module 04): lags beyond the record wrap around. Zero-pad both records to length \(\ge N_x + N_y - 1\) to obtain the linear correlation. For kernel \(W = 80\) and search ±24 lags, a direct sum costs 80 × 49 ≈ 3.9k multiply–adds per estimate; FFT correlation wins only for long windows or many lags (cost \(\propto N\log N\)).

Demo 11.1

NCC delay estimate on two speckle segments

SimulationSeeded speckle synthesized from a Gaussian pulse spectrum (FFT, JS), shifted by \(D\), partly decorrelated and noised; NCC (11.3) and parabolic fit computed in the page.
x(t) reference (kernel window shaded)y(t) delayed + noise + decorrelation
NCC samples ρ̂(k)parabolic fittrue shiftestimate
Try: \(f_s\) = 40 MHz, kernel 2 µs (80 samples). Set SNR to 0 dB and \(\rho\) to 0.6, then reseed a few times: the global maximum sometimes jumps to a neighbouring peak one RF period (\(1/f_0\)) away. Narrow the bandwidth to 20 % and the side peaks rise toward the main peak, so these jumps get more frequent.
Ultrasound link

Two-way geometry: \(\Delta z = cD/2\), so 1 ns of delay is 0.77 µm of axial displacement at \(c\) = 1540 m/s, and one sample at \(f_s\) = 40 MHz (25 ns) is 19.3 µm. ARFI displacements are 1–20 µm, which means 0.05–1 sample. Every practical tracker is therefore a subsample estimator, and its bias (Section 11.2) matters as much as its jitter (Section 11.3).

Pitfall

A high peak NCC does not guarantee a correct delay. With narrowband RF the side peaks at \(\pm 1/f_0\) can be 0.8 or higher. NCC tells you the two windows are similar at the chosen lag, not that you chose the right lobe. Use envelope or IQ-magnitude correlation for a coarse search and RF or phase for refinement, or regularize with neighbours.

11.2Subsample estimation and its bias

The NCC is only known at integer lags. Let \(k_0\) be the integer peak and \(R_-\), \(R_0\), \(R_+\) the values at \(k_0 - 1\), \(k_0\), \(k_0 + 1\). Each interpolator fits a model to these three points (or to more) and returns \(k_0 + \delta\).

EstimatorFractional offset \(\delta\)Model assumed for the peak
Parabolic\(\delta = \dfrac{R_- - R_+}{2(R_- - 2R_0 + R_+)}\)2nd-order Taylor expansion; exact only for a quadratic peak
Cosine (Céspedes et al. 1995)\(\omega = \arccos\dfrac{R_- + R_+}{2R_0}\), \(\delta = \dfrac{1}{\omega}\arctan\dfrac{R_+ - R_-}{2R_0\sin\omega}\)\(R(k) \propto \cos(\omega(k - \delta))\): exact for a narrowband RF correlation
Gaussianparabolic fit applied to \(\ln R\)Gaussian peak: exact for the magnitude \(\lvert C(\tau)\rvert\) of the complex correlation of a Gaussian spectrum, close for envelope correlation, poor for RF
Upsample ×U (sinc)band-limited interpolation of \(R\) on a \(1/U\) grid, then a parabolic fit on that fine grid\(R\) band-limited (true when \(f_s\) exceeds twice the highest RF frequency); the residual parabolic bias falls as \(U^{-3}\) in original samples (0.012 → 0.0014 → 0.00018 for U = 1, 2, 4), until the finite interpolation kernel sets a floor (\(3.5\times10^{-4}\) sample here with a 12-lag Hann-tapered sinc)
Phase zero-crossing\(\delta = -\angle C(k_0)/(2\pi f_0 T_s)\), with \(C\) the complex correlation of the analytic signals (Module 08)\(C(\tau) = R_{\mathrm{env}}(\tau)\, e^{j2\pi f_0\tau}\) with \(R_{\mathrm{env}}\) real when the spectrum is symmetric about \(f_0\): the phase is exactly linear in \(\tau\), so there is no bias if \(f_0\) equals the spectral centroid

Why the bias is zero at 0 and 0.5 samples

\[ \begin{gathered} R(k) = R_{ss}\bigl((k - d)\,T_s\bigr), \quad R_{ss} \text{ even} \\ d = 0:\ R_- = R_+ \;\Rightarrow\; \delta = 0 \qquad d = \tfrac{1}{2}:\ R(0) = R(1) \;\Rightarrow\; \text{any symmetric fit returns } \tfrac{1}{2} \end{gathered} \]
(11.8)
\(R(k)\)
correlation function sampled at integer lag \(k\)
\(d\)
true fractional delay in samples, \(0 \le d < 1\)
\(R_{ss}\)
expected (noise-free, long-window) correlation of the speckle, even in its argument
\(T_s\)
sampling interval, s
\(R_-,\ R_+\)
correlation values at the lags \(-1\) and \(+1\) around the integer peak
\(\delta\)
fractional offset returned by the subsample fit, samples
\[ \mathrm{bias}(d) = -\,\mathrm{bias}(1 - d) \;\Rightarrow\; \mathrm{bias}(d) \text{ is odd about } d = \tfrac{1}{2}, \text{ with zeros at } d \in \bigl\{0, \tfrac{1}{2}, 1\bigr\} \]
(11.9)
\(\mathrm{bias}(d)\)
subsample estimator bias at true fractional delay \(d\), samples
\(d\)
true fractional delay in samples, as in (11.8)
The extremes sit in between, near \(d \approx 0.3\) and 0.7 for the RF fits here. The parabolic and cosine fits underestimate \(\lvert d - k_0\rvert\), so the estimates cluster toward integer lags. That is pixel locking (in PIV and speckle-tracking terms, "peak locking"), and it shows up as a sawtooth-like texture in displacement maps. The Gaussian fit applied to RF errs the other way (toward half-samples), with twice the parabolic magnitude.

The parabolic bias depends on the peak curvature relative to the sampling interval, so on \(f_0/f_s\). With \(f_0\) = 5 MHz and \(f_s\) = 40 MHz, the RF correlation completes a cycle in 8 samples and the parabola underestimates \(\lvert\delta\rvert\) by about 0.012 samples (0.3 ns, 0.23 µm) at worst. That is small, but it is systematic: averaging does not remove it, and it does not cancel when you take spatial derivatives for strain.

Demo 11.2

Bias curves of five subsample estimators

Approx. modelEach estimator is run on the expected RF correlation of Gaussian-spectrum speckle (noise-free, infinite window), sampled at \(k - d\) for 41 values of \(d\) (JS). Real kernels scatter around these curves.
paraboliccosineGaussianupsample ×8 (sinc)phase (complex corr.)
Try: drop \(f_s\) to 20 MHz at \(f_0\) = 5 MHz: the parabolic bias grows from 0.012 to 0.052 sample (4.4×, about \((f_0/f_s)^2\)), 2.6 ns. Raise \(f_0\) to 7.5 MHz at 20 MHz (\(f_s/f_0\) = 2.7) and the integer maximum lands on the wrong RF lobe for part of the \(d\) range: every interpolator then fails by whole periods, which the readout flags. Then set the \(f_0\) error to +10 % (the assumed frequency is 10 % above the echo's centre frequency, as after an attenuation downshift): the phase estimator, unbiased until now, picks up a bias \((k_0 - d)\,e/(1+e)\) proportional to the distance from the integer lag, up to 0.045 sample, and its curve jumps at \(d = \tfrac{1}{2}\) instead of crossing zero.
Intuition

An RF correlation peak looks like a cosine under a Gaussian envelope. A parabola fitted through three points of a cosine is too flat at the top, so it pulls the estimate toward the centre sample. The cosine fit has the right shape and is nearly exact; residual bias comes only from the envelope. The phase estimator skips the shape question: it reads the delay from the angle of the complex correlation, which is linear in \(\tau\).

Pitfall

The phase estimator converts angle to delay with an assumed \(f_0\). Frequency-dependent attenuation lowers the received centre frequency with depth (Module 08: at 0.5 dB/cm/MHz and 60 % bandwidth, 5 MHz drops by 0.75 MHz at 4 cm), so a fixed \(f_0\) gives a depth-dependent scale error. Loupas' 2D autocorrelator (Section 11.4) estimates the local mean RF frequency from the same data and removes this error.

Ultrasound link

Pinton, Dahl & Trahey (2006) compared normalized cross-correlation and phase-shift estimators for ARFI-scale displacements (below half a wavelength). Loupas' 2D autocorrelator generally matched NCC accuracy at a fraction of the computational cost. Upsampled or sinc-interpolated NCC is the bias-free reference you would use for a phase-screen ground truth.

11.3Variance limits: the Walker–Trahey bound and the threshold effect

Bias is systematic; jitter is random. For two jointly Gaussian records with coherence \(\gamma(f)\), the Cramér–Rao bound on any unbiased delay estimate from a window of duration \(T\) is (Carter 1987, building on Knapp & Carter 1976):

\[ \sigma_D^2 \;\ge\; \Bigl[\, 2T \int_0^\infty (2\pi f)^2\, \frac{\lvert\gamma(f)\rvert^2}{1 - \lvert\gamma(f)\rvert^2}\, df \,\Bigr]^{-1} \]
(11.10)
\(\sigma_D\)
standard deviation of the delay estimate, s
\(\gamma(f)\)
complex coherence between the two records at frequency \(f\)
\(f\)
frequency, Hz
\(T\)
observation window, s
The Fisher information grows with window length \(T\), with the second moment of frequency (high frequencies give sharp peaks), and with \(\lvert\gamma\rvert^2/(1 - \lvert\gamma\rvert^2)\), which diverges as the records become identical.

Walker & Trahey (1995) applied this to speckle: a rectangular band of centre \(f_0\) and fractional width \(B\), with the coherence set by both the electronic SNR and the speckle correlation \(\rho\):

\[ \begin{gathered} \lvert\gamma\rvert = \rho \cdot \frac{\mathrm{SNR}_p}{1 + \mathrm{SNR}_p} \quad\Rightarrow\quad \frac{1}{\lvert\gamma\rvert^2} - 1 = \frac{1}{\rho^2}\Bigl(1 + \frac{1}{\mathrm{SNR}^2}\Bigr)^2 - 1 \\ (\mathrm{SNR} = \text{amplitude ratio}, \ \mathrm{SNR}_p = \mathrm{SNR}^2) \end{gathered} \]
(11.11)
\(\lvert\gamma\rvert\)
coherence magnitude over the band
\(\rho\)
speckle correlation, as in (11.6)
\(\mathrm{SNR}\)
electronic signal-to-noise amplitude ratio
\(\mathrm{SNR}_p\)
signal-to-noise power ratio, \(\mathrm{SNR}^2\)
\[ \int_{f_0(1 - B/2)}^{f_0(1 + B/2)} (2\pi f)^2\, df = 4\pi^2 f_0^3 \Bigl(B + \frac{B^3}{12}\Bigr) = \frac{\pi^2}{3}\, f_0^3\, \bigl(B^3 + 12B\bigr) \]
(11.12)
\(f\)
frequency, Hz
\(f_0\)
centre frequency, Hz
\(B\)
fractional bandwidth (band width divided by \(f_0\)), dimensionless
\((2\pi f)^2\)
squared angular frequency weight from (11.10), rad²/s²
\[ \sigma_D \;\ge\; \sqrt{ \frac{3}{2\pi^2 f_0^3\, T\, \bigl(B^3 + 12B\bigr)} \;\Bigl[ \frac{1}{\rho^2}\Bigl(1 + \frac{1}{\mathrm{SNR}^2}\Bigr)^2 - 1 \Bigr] } \]
(11.13)
\(\sigma_D\)
lower bound on the delay standard deviation, s
\(f_0\)
centre frequency, Hz
\(B\)
fractional bandwidth, dimensionless
\(T\)
window duration, s
\(\rho\)
speckle correlation, 0 to 1
\(\mathrm{SNR}\)
electronic signal-to-noise amplitude ratio
This is the Walker–Trahey bound: the Carter form evaluated with flat coherence across the band, derived by (11.10)–(11.12). \(f_0\) in Hz and \(T\) in s give \(\sigma_D\) in s. Example: \(f_0\) = 5 MHz, \(B\) = 0.6, \(T\) = 2 µs, SNR = 30 dB, \(\rho\) = 1 gives \(\sigma_D \ge\) 0.41 ns (0.31 µm); \(\rho\) = 0.99 gives 1.35 ns.

The bound shows four scalings: \(\sigma \propto f_0^{-3/2}\), \(\propto T^{-1/2}\), \(\approx \propto B^{-1/2}\) (for \(B \lesssim 1\)), and at high SNR a floor \(\sigma \propto \sqrt{1/\rho^2 - 1} \approx \sqrt{2(1-\rho)}\) set by decorrelation. Once \(\rho < 1\), raising the SNR stops helping beyond the point where \(1/\mathrm{SNR}_p \approx 1 - \rho\).

The threshold effect

The CRLB is a local bound: it describes small errors around the true peak. The RF correlation has side peaks at \(\tau \approx \pm m/f_0\). When noise or decorrelation makes a side peak exceed the main peak, the error jumps by a whole RF period (200 ns, or 154 µm in displacement at 5 MHz). The error distribution becomes a mixture: CRLB-like jitter plus outliers spaced \(1/f_0\) apart, ending in a uniform distribution over the search range. Below a threshold SNR (or \(\rho\), or \(T\cdot B\)) the variance rises far above the CRLB. Global bounds that account for these "ambiguity" errors, the Ziv–Zakai and Barankin bounds (Weiss & Weinstein 1983), reproduce this knee. Narrow bandwidth, short kernels and large search ranges all move the threshold to higher SNR.

Demo 11.3

Monte Carlo jitter vs SNR against the CRLB

Simulation300 seeded Monte Carlo trials per SNR: speckle pairs synthesized in the frequency domain, NCC + parabola and complex-correlation phase estimators (JS); CRLB from the Carter integral (11.10) for the same spectrum, numerically integrated.
NCC + parabolic (MC)phase of complex corr. (MC)CRLB, Carter integral for this spectrumWalker–Trahey closed form
Try: \(f_0\) = 5 MHz, \(f_s\) = 40 MHz, true shift uniform in [0, 1) sample per trial, SNR from −10 to 40 dB in 5 dB steps. From 15 to 30 dB both estimators follow the bound to within about 25 %. At 10 dB about 5 % of trials hop and \(\sigma\) is roughly 18× the CRLB; at 5 dB more than a third hop. That knee is the threshold. Shorten \(T\) to 0.5 µs (2.5 periods): hops reach about a quarter of trials at 15 dB and 1–3 % remain even at 40 dB, because a window that short looks almost sinusoidal and neighbouring RF cycles match it nearly as well. Narrow the bandwidth to 20 % (side peaks at 0.93) and about 9 % of trials still hop at 30 dB. Set \(\rho\) = 0.99 and both curves flatten into the decorrelation floor (the Carter bound includes \(\rho\), so the ratio stays near 1.1–1.5). At very high SNR with \(\rho\) = 1 the simulation sits above the bound: parabolic bias and window-edge decorrelation (the kernels overlap only on \(W - \lvert D\rvert\) samples) add a floor of about 0.3 ns. The Walker–Trahey closed form assumes a rectangular band: for this Gaussian spectrum it matches the Carter integral at \(\rho\) = 1 and is about 20 % higher at \(\rho\) = 0.99.
Ultrasound link

The bound covers all three targets. Phase-screen estimation uses neighbouring elements with \(\rho\) set by the van Cittert–Zernike coherence (Module 07). For diffuse scatter at the transmit focus, adjacent elements have \(\rho = 1 - 1/N\) for a uniform \(N\)-element transmit (0.984 for \(N\) = 64), lower away from focus or through an aberrator, so \(\rho \approx\) 0.9–0.98 and the floor term dominates. Sound-speed methods average many such delays. ULM localizes isolated bubbles with high SNR and \(\rho \approx 1\), which is why it can beat \(\lambda\) by a factor of 10.

Pitfall

Do not report "\(\sigma\) = CRLB" from simulations that never leave the high-SNR regime. Real data contain a mix of good and threshold-region kernels. Report the outlier fraction separately, or use a robust spread (median absolute deviation) together with the plain standard deviation.

11.4Phase-based estimators on IQ: Kasai and Loupas

After IQ demodulation (Module 08), a scatterer moving axially at \(v\) rotates the phase of the slow-time sample at a fixed depth by a constant step per pulse:

\[ z[m] = A\, e^{j\varphi_0}\, e^{\,j m \Delta\varphi} + n[m], \qquad \Delta\varphi = 2\pi f_0 \cdot \frac{2v}{c} \cdot T_{\mathrm{PRF}} = \frac{4\pi f_0 v}{c \cdot \mathrm{PRF}} \]
(11.14)
\(z[m]\)
slow-time IQ sample of pulse \(m\) at one depth gate
\(A,\ \varphi_0\)
amplitude and initial phase of the echo, rad
\(m\)
pulse (slow-time) index
\(n[m]\)
additive noise
\(\Delta\varphi\)
phase step per pulse, rad
\(f_0\)
transmit centre frequency, Hz
\(v\)
axial velocity, m/s
\(c\)
sound speed, 1540 m/s
\(\mathrm{PRF}\)
pulse repetition frequency, Hz; \(T_{\mathrm{PRF}} = 1/\mathrm{PRF}\)
\(j\)
imaginary unit
\[ \hat{R}(1) = \sum_g \sum_{m=0}^{M-2} z_g[m+1]\, z_g^*[m], \qquad \hat{v} = \frac{c \cdot \mathrm{PRF}}{4\pi f_0} \cdot \angle \hat{R}(1) \]
(11.15)
\(\hat{R}(1)\)
lag-one slow-time autocorrelation
\(z_g[m]\)
slow-time samples of the \(g\)th gate, as in (11.14)
\((\cdot)^*\)
complex conjugate
\(M\)
ensemble length (pulses)
\(g\)
index of the averaged gates (spatial or channel)
\(\hat{v}\)
axial velocity estimate, m/s
\(c,\ f_0,\ \mathrm{PRF}\)
sound speed, centre frequency and pulse repetition frequency, as in (11.14)
\(\angle\)
phase angle, rad, wrapped to \((-\pi, \pi]\)
This is the Kasai lag-one autocorrelation estimator (1985). \(\angle R(1)/(2\pi T_{\mathrm{PRF}})\) is the mean Doppler frequency (the centroid of the slow-time power spectrum, exactly for a symmetric spectrum). It costs one complex multiply–add per sample. It is not the ML frequency estimator: for a stationary tone at high SNR the sum telescopes, so it effectively uses only the first and last pulses, and its variance exceeds the CRLB by 1.8× at \(M\) = 8 and 14× at \(M\) = 64 (m11_kasai_doppler.m, 20 dB; the high-SNR limit is \(M(M+1)/(6(M-1))\), 1.7× and 11×).

\(\angle\hat{R}(1)\) is only defined modulo \(2\pi\), so velocities alias when \(\lvert\Delta\varphi\rvert > \pi\):

\[ v_{\mathrm{Nyq}} = \frac{c \cdot \mathrm{PRF}}{4 f_0} \qquad \text{(displacement per pulse limited to } \pm\lambda/4) \]
(11.16)
\(v_{\mathrm{Nyq}}\)
Nyquist (aliasing) velocity, m/s
\(c,\ f_0,\ \mathrm{PRF}\)
as in (11.14)
\(\lambda\)
wavelength \(c/f_0\), m
\(f_0\) = 5 MHz, PRF = 5 kHz: \(v_{\mathrm{Nyq}} = 1540 \cdot 5000 / (2\times10^7)\) = 0.385 m/s. The same \(\pm\lambda/4\) wrap limits phase-based displacement tracking between any two acquisitions, so ARFI uses PRFs of 5–10 kHz.

Loupas 2D autocorrelation (1995) also takes the lag-one correlation along depth, \(R\)(fast-time lag 1, slow-time lag 0), to estimate the local mean RF frequency \(\bar{f}\). It then uses \(\hat{v} = (c\cdot\mathrm{PRF}/4\pi)\cdot\angle R(0,1)/\bar{f}\). This removes the \(f_0\)-mismatch bias shown in Demo 11.2 and is a common choice for ARFI/SWE displacement tracking.

Speed–accuracy tradeoff. A higher PRF raises \(v_{\mathrm{Nyq}}\) but makes \(\Delta\varphi\) per pulse smaller, so for fixed phase noise the velocity jitter \(\sigma_v \propto \mathrm{PRF}\cdot\sigma_\varphi\) grows. A longer ensemble \(M\) lowers \(\sigma_\varphi\), as \(1/(M-1)\) for a stationary tone at high SNR and roughly as \(M^{-1/2}\) for a random Doppler signal with a finite spectral width (blood, decorrelating speckle), but it costs frame rate and admits more decorrelation.

PropertyTime-domain (NCC + interpolation)Phase-domain (Kasai / Loupas on IQ)
Biasinterpolator bias, ~0.01 sample (parabolic), ≪ with cosine or upsamplingnone if \(f_0\) (Kasai) or \(\bar{f}\) (Loupas) is right; scale error \(\propto f_0\) mismatch
Variancenear CRLB above thresholdnear CRLB above threshold for small shifts; similar to NCC
Aliasing / ambiguityfalse peaks at \(\pm 1/f_0\) (peak hopping); range limited by search window, not by \(\lambda/4\)hard wrap at \(\pm\lambda/4\) per lag; needs unwrapping or a coarse NCC stage
Large displacementsyes (search \(\pm K\) lags)no (\(\lvert D\rvert < 1/(2f_0)\))
Cost per estimate\(W\cdot(2K+1)\) MACs (≈4k for \(W\) = 80, \(K\) = 24) or FFT≈ \(W\) complex MACs (≈80)
Typical usespeckle tracking, large strain, phase screens with large delays, ULM patch matchingcolour Doppler, ARFI, SWE, small-displacement tracking
Demo 11.4

Kasai velocity estimator and aliasing

SimulationSeeded slow-time ensembles (constant phase step per pulse plus complex white noise, 8 gates) and the estimator (11.15), computed in the page (JS).
theory: wrapped vKasai estimates (one ensemble each)current v
Try: \(f_0\) = 5 MHz, 8 axial gates. Push \(v\) past \(v_{\mathrm{Nyq}}\): in the lower plot the phasor step exceeds 180°, so the estimator reads it as a smaller step the other way, and the estimate jumps to the opposite sign. Lower the PRF to 2 kHz and \(v_{\mathrm{Nyq}}\) falls to 0.15 m/s. Lower the SNR to 0 dB with \(M\) = 3 and the scatter grows near \(\pm v_{\mathrm{Nyq}}\), where wrap errors flip the sign.
Ultrasound link

In ULM, bubbles are separated from tissue by SVD clutter filtering (Module 06), then localized frame by frame. Phase-based velocity is rarely used for the bubbles themselves, because each bubble is a point target, not stationary speckle. Tissue motion accumulated over a ULM acquisition (breathing, pulsation, seconds of data) spans many wavelengths, far beyond the \(\pm\lambda/4\) wrap, so motion correction usually registers images (NCC or phase correlation of B-mode blocks) rather than tracking phase.

11.5Decorrelation sources and kernel length

\(\rho\) in the bound collects every way the second window differs from a shifted copy of the first:

SourceMechanismScaling / rule of thumb
Electronic noiseindependent noise in each record\(\lvert\gamma\rvert = \mathrm{SNR}_p/(1+\mathrm{SNR}_p)\); enters as \((1 + 1/\mathrm{SNR}^2)^2\)
Axial strain \(\varepsilon\)window is compressed: local shift varies as \(\varepsilon t\) across the kernel\(\rho(\varepsilon, T) \approx \frac{1}{T}\int_{-T/2}^{T/2} R_n(\varepsilon t)\, dt\); collapses once \(\varepsilon T \gtrsim 1/(2f_0)\)
Lateral / elevational motionscatterers leave the PSF\(\rho \approx\) normalized autocorrelation of the lateral PSF; Gaussian PSF of amplitude FWHM \(w\): \(\rho = \exp[-2\ln 2\,(\Delta x/w)^2]\), so \(\rho\) = 0.92 at \(w/4\), 0.71 at \(w/2\), 0.25 at \(w\)
Shear, rotationdifferent parts of the kernel move differentlylike strain, but in the lateral direction
Bandwidthnarrow band: broad envelope, but high side peakslower threshold SNR at large \(B\); \(f_0^3 B\) scaling of the Fisher information
Element-to-element (aberration)different paths through tissue and different spatial positionsvan Cittert–Zernike: \(\rho(\Delta)\) is the autocorrelation of the transmit apodization

Kernel length is a three-way tradeoff. A longer \(T\) reduces jitter as \(T^{-1/2}\), but (i) axial resolution of the displacement or strain map worsens to about \(cT/2\), and (ii) under strain the decorrelation term grows with \(T\). The combined \(\sigma(T)\) has a minimum: the kernel analogue of the elastographic strain filter (Varghese & Ophir 1997).

Demo 11.5

Kernel length vs strain: finding the optimum window

Approx. modelWalker–Trahey bound (11.13) with \(\rho(\varepsilon, T)\) = window average of the expected Gaussian-pulse RF correlation under uniform strain, no companding; ignores the threshold effect.
ε = 0 (∝ T−1/2)with strain ερ(ε, T) (right axis)optimum
Try: at \(\varepsilon\) = 1 % and 20 dB SNR the optimum is about 1.5 µs (\(\sigma \approx\) 2.1 ns). At 3 % it shrinks to about 0.5 µs (\(\sigma \approx\) 3.6 ns). Raise the SNR to 40 dB and the optimum shortens further, because the noise term no longer needs a long window to average down. The bound ignores the threshold effect, so kernels shorter than about two RF periods (shaded) are not usable in practice. This is why large-strain elastography uses companding (pre-stretching one record) or multi-step tracking.
Intuition

Under strain, a long kernel compares an accordion with a stretched accordion: the centres line up while the ends drift apart by \(\varepsilon T/2\). The bound balances two terms in its bracket: the noise term \(1/\mathrm{SNR}_p\) and the decorrelation term \(1 - \rho^2\), which grows as \((\varepsilon T f_0)^2\). The optimum sits where they are comparable, so it lies well before the ends drift by a full period. At 20 dB SNR it falls near \(\varepsilon T \approx 0.075/f_0\) (end drift ≈ 4 % of a period); at 30 dB, near \(0.025/f_0\).

11.6Ultrasound applications

(a) Near-field phase-screen estimation (Flax & O'Donnell 1988)

Model the aberrator as a thin time-delay screen \(\tau_i\) at the aperture. Flax & O'Donnell cross-correlated signals from adjacent elements, from point reflectors or diffuse speckle, to estimate the differences \(\Delta\tau_i = \tau_{i+1} - \tau_i\). They integrated these into a profile and applied it as a correction, iterating because the transmit focus improves too. Adjacent elements are used because their speckle is highly correlated (van Cittert–Zernike, Module 07), while distant elements are not.

SchematicGeometry of a near-field phase screen, not to scale (depth compressed, screen thickness exaggerated). Typical values: 0.3 mm pitch (≈ \(\lambda\) at 5 MHz), screen 20–80 ns rms with 3–6 mm correlation length (a fat–muscle body wall), target 20–60 mm deep. The right panel is an illustrative estimate, not a computed run (Demo 11.6 computes one).
Thin delay screen at the aperture aberrating layer: τi ∝ thickness × (1/clayer − 1/c) array, pitch 0.3 mm NCC(i, i+1) → Δτ̂i point target or speckle region, z ≈ 20–60 mm wavefront after the screen: delayed by τi per element τ (ns) vs element n — true τn‒ ‒ Σ Δτ̂, drifts as √n σΔ
\[ \hat{\tau}_n = \sum_{i=0}^{n-1} \Delta\hat{\tau}_i \quad\Rightarrow\quad \mathrm{Var}\bigl(\hat{\tau}_n - \tau_n\bigr) = n\, \sigma_\Delta^2 \qquad \text{(independent errors: a random walk)} \]
(11.17)
\(\hat{\tau}_n\)
accumulated delay estimate at element \(n\), s
\(\tau_n\)
true delay at element \(n\), s
\(\Delta\hat{\tau}_i\)
estimated delay between neighbouring elements \(i\) and \(i+1\), s
\(\sigma_\Delta\)
standard deviation of each neighbour estimate, s
\(n,\ i\)
element indices
\(\mathrm{Var}\)
variance
With \(\sigma_\Delta\) = 2 ns per pair, the error at element 64 is about 16 ns (≈ 29° at 5 MHz). A single peak hop adds a \(1/f_0\) step to every element after it. Cures: correlate each element against a beamsum reference (errors no longer accumulate), least-squares integration of multiple lags (\(\Delta\tau\) from \(i\) to \(i+1\), \(i+2\), …), and NCC-weighted integration.

The constant and linear terms of \(\tau_i\) are not observable from diffuse speckle (they shift or steer the image without defocusing it), so compare profiles after removing the mean and linear trend.

Demo 11.6

Element-to-element delay estimation across a 64-element aperture

SimulationSeeded Gaussian-correlated random screen; per-element speckle with adjacent correlation \(\rho\) and noise synthesized in the frequency domain; adjacent NCC + parabola and cumulative sum (JS).
true screen τicumulative sum of adjacent NCC delays
error τ̂n − τn (anchored at element 0)±√n σΔ random-walk envelope (CRLB)
Try: \(f_0\) = 5 MHz, \(B\) = 60 %, \(f_s\) = 40 MHz, 2 µs kernel, pitch 0.3 mm. The error curve wanders like a random walk inside the \(\pm\sqrt{n}\) envelope. Lower the SNR to 0 dB or \(\rho\) to 0.7: some adjacent estimates hop a full period (200 ns), and every later element inherits the step. Raise the screen rms with a short correlation length and adjacent delays approach the \(\pm K\) search limit.

(b) Speckle tracking and elastography

Quasi-static elastography estimates displacement \(u(z)\) with NCC or phase tracking and differentiates it to get strain, \(\varepsilon \approx \Delta u/\Delta z\). Differentiation amplifies jitter: for two independent estimates spaced \(\Delta z\) apart, \(\sigma_\varepsilon \approx \sqrt{2}\,\sigma_u/\Delta z\). With \(\sigma_u\) = 0.3 µm and \(\Delta z\) = 1 mm, \(\sigma_\varepsilon \approx\) 0.04 % strain. ARFI and SWE track 1–20 µm displacements at PRFs of 5–10 kHz with Loupas or Kasai. For SWE on the uML15-3F class of probes, the shear-wave arrival time is itself a time-delay estimate between lateral positions, and the same bias/variance machinery applies to cross-correlating particle-velocity traces.

(c) Sound-speed estimation from delays and phase shifts

Computed ultrasound tomography in echo mode (CUTE; Jaeger et al. 2015) beamforms the same region with several transmit steering angles under an assumed \(c\). If the true \(c\) differs locally, echoes from a given point arrive with angle-dependent delay errors. The phase shift between complex (IQ) images from neighbouring angles, estimated with a Kasai-type lag-one product across angles, is a line integral of the slowness error along the two paths. A linear inverse problem (Batch 2) then recovers the local \(c\). The phase estimator's \(\pm\lambda/4\) wrap and its noise set the usable angle spacing; at 5 MHz, 1 ns corresponds to 1.8° of phase.

(d) ULM bubble tracking and pairing

After localization, ULM links the bubble positions in frame \(k\) to those in frame \(k+1\). Two families of methods: (i) NCC of small image patches around each detection, the 2D version of 11.1, and (ii) global assignment that minimises total travel distance, solved with the Hungarian (Kuhn–Munkres) algorithm or a Kalman-filter variant. At 1 kHz and 10 mm/s the displacement per frame is 10 µm (≈ \(\lambda/30\) at 5 MHz). The search radius must exceed the maximum flow displacement but stay small enough to avoid mispairing in dense bubble populations. That is the same threshold effect as in 11.3, now in assignment space. See Target A (ULM).

Ultrasound link

A 128-element, 0.3 mm pitch array at 5 MHz (\(\lambda\) = 0.31 mm) samples the aperture at about \(\lambda\). Flax & O'Donnell-style integration over 128 elements gives \(\sqrt{128} \approx 11\) times the adjacent-pair jitter at the far edge. Distortion-matrix and reflection-matrix methods (Module 33) instead estimate the screen from the whole set of element pairs at once.

Pitfall

Adjacent-element correlation measures the difference of the screen. A large-scale sound-speed error (a smooth quadratic delay profile) has tiny adjacent differences that are easily dominated by bias. Such focusing errors are better estimated with a global metric (coherence factor, Module 07) or with tomographic sound-speed estimation.

11.7Visual summary

Ten pictures that carry the module, each computed from the same correlation, bound and Monte Carlo models as the demos above. Each card is a question you should answer in five seconds, the answer with its number, and the figure to keep in mind. The slider sets the fractional bandwidth used by the correlation-shape cards: narrow bands give broad envelopes but tall false peaks.

11.8How experts think about it

Every delay estimator is a peak finder on a correlation function shaped like a cosine at \(f_0\) under an envelope of width ~\(1/B\). Experts separate three error types and size each before trusting a number: bias (interpolator shape mismatch, wrong \(f_0\)), jitter (the local CRLB, set by \(f_0\), \(B\), \(T\), SNR and \(\rho\)), and false peaks (whole RF periods, which no bound describes). Then they convert to µm with the two-way factor and ask which of the three dominates.

The correlation peak is the maximum-likelihood delay in white Gaussian noise.
matched filterGCC weight \(\lvert\gamma\rvert^2/(\lvert G_{xy}\rvert(1-\lvert\gamma\rvert^2))\)

With a known template \(s\) and white Gaussian noise of variance \(\sigma^2\), \(\ln p(y; D)\) \(= -\frac{1}{2\sigma^2}\sum (y[n] - s[n-D])^2\) + const (11.4). The energy term \(\sum s^2[n-D]\) does not depend on \(D\), so maximizing the likelihood means maximizing \(\sum y[n]\,s[n-D]\): the cross-correlation. With both records noisy (speckle tracking), Knapp & Carter's ML processor weights the cross-spectrum first; for white noise and in-band SNR ≫ 1 that weight is flat across the band, so correlation of band-passed RF is near ML.

SimulationOne seeded noisy realization of a 5 MHz, 60 % pulse at 40 MHz; both curves computed over trial delays (JS).
Coloured noise or clutter (pre-whiten or filter first); unequal spectra; low SNR, where the ML estimate is dominated by false peaks and is no longer near-unbiased; signal that changes shape between records (strain, decorrelation), which the shift model ignores.
One nanosecond of delay is 0.77 µm of axial displacement; one sample at 40 MHz is 19.25 µm.
\(\Delta z = cD/2\)1 ns ↔ 0.77 µm1 ns ↔ 1.8° at 5 MHz

Pulse-echo delay covers the path twice, so \(\Delta z = cD/2 = 1540 \times 10^{-9}/2\) = 0.77 µm. At \(f_s\) = 40 MHz, \(T_s\) = 25 ns ↔ 19.25 µm, so ARFI displacements of 1–20 µm are 0.05–1 sample: every tracker is a subsample estimator. Phase: 360° × \(f_0\) × 1 ns = 1.8° at 5 MHz.

Exact formula
Assumes the delay is a round trip at one sound speed. Element-to-element receive delays in phase-screen work are one-way (1 ns ↔ 1.54 µm of path); transmit-angle delays in sound-speed tomography mix both.
Parabolic interpolation is biased toward the integer sample; the cosine fit is nearly unbiased on RF; any symmetric fit is exact at 0 and ½ sample.
parabolic 0.012 smp at \(f_s/f_0\) = 80.052 at \(f_s/f_0\) = 4cosine 0.0014

\(R(k) = R_{ss}((k-d)T_s)\) with \(R_{ss}\) even, so \(d = 0\) gives \(R_- = R_+\) and \(d = \tfrac{1}{2}\) gives \(R(0) = R(1)\): symmetric fits return 0 and ½ exactly, and \(\mathrm{bias}(d) = -\mathrm{bias}(1-d)\) (11.9). In between, a parabola through three points of a cosine is too flat at the top and pulls toward \(k_0\): worst 0.012 sample (0.30 ns) at 5 MHz, 40 MHz, rising about as \((f_0/f_s)^2\). The cosine model matches the carrier and leaves only envelope error.

Approx. modelFits run on the expected correlation of Gaussian-spectrum speckle (5 MHz, 60 %), noise-free, infinite window.
The Gaussian fit on RF errs the other way (0.025 sample, toward half-samples). The curves are expectations: one noise-free kernel's error varies by ±0.011 sample with its local spectrum (m11_ncc_subsample_bias.m), so you cannot subtract the mean curve. Below \(f_s/f_0 \approx 3\) the integer peak can land on the wrong lobe.
For sub-nanosecond precision, upsample RF to at least ~\(8f_0\), or work on the phase of complex baseband.
\(f_s/f_0\) = 4: 2.6 ns bias\(f_s/f_0\) = 8: 0.30 nsupsample ×U: bias \(\propto U^{-3}\)

Parabolic bias scales with the peak curvature per sample, \((2\pi f_0/f_s)^2\). At 5 MHz, \(f_s\) = 20 MHz leaves 0.052 sample × 50 ns = 2.6 ns; \(f_s\) = 40 MHz leaves 0.30 ns. Band-limited (sinc) upsampling of the correlation by \(U\) before the parabola cuts the residual as \(U^{-3}\) in original samples (0.012, 0.0014, 0.00018 for U = 1, 2, 4). The phase of the complex correlation needs no shape model at all, so a modest IQ rate suffices as long as the residual after the integer search stays inside \(\pm 1/(2f_0)\), i.e. an IQ sample period shorter than \(1/f_0\).

Approx. modelMaximum parabolic bias on the expected Gaussian-spectrum correlation, swept over \(f_s/f_0\).
Upsampling cannot create information: it fixes bias, not jitter. A short interpolation kernel sets a floor (\(3.5\times10^{-4}\) sample for a 12-lag tapered sinc). The phase route needs the right \(f_0\).
Phase-based estimators on IQ are cheap and nearly unbiased, but they wrap at \(\pm\lambda/4\) of displacement.
\(\lvert D\rvert < 1/(2f_0)\) = 100 ns at 5 MHz±77 µm~\(W\) vs \(W(2K+1)\) MACs

The complex correlation of analytic signals is \(C(\tau) = R_{\mathrm{env}}(\tau)\,e^{j2\pi f_0\tau}\), so \(\hat{D} = -\angle C/(2\pi f_0)\) with no interpolation model. \(\angle C\) is known only modulo \(2\pi\), so the unambiguous range is ±half a period: ±100 ns of delay, \(\pm\lambda/4\) = ±77 µm of displacement at 5 MHz. Cost: one complex multiply–add per sample (≈80) versus 80 × 49 ≈ 3.9k for a ±24-lag NCC search.

Exact formula
An \(f_0\) error \(e\) gives a bias \((k_0 - d)\,e/(1+e)\), up to 0.024 sample for 5 %; attenuation downshift makes this depth dependent (Loupas estimates the local \(\bar{f}\)). Large motion needs a coarse NCC stage or unwrapping.
Jitter follows the Walker–Trahey bound: \(\sigma\) falls as \(f_0^{-3/2}\), \(T^{-1/2}\) and roughly \(B^{-1/2}\).
0.41 ns at 5 MHz, 60 %, 2 µs, 30 dB1.35 ns at \(\rho\) = 0.99

Exact form used on this page, (11.13) (rectangular band of fractional width \(B\), SNR as amplitude ratio):

Exact formula
\[ \sigma_D \ge \sqrt{ \frac{3}{2\pi^2 f_0^3 T (B^3 + 12B)} \Bigl[ \frac{(1 + 1/\mathrm{SNR}^2)^2}{\rho^2} - 1 \Bigr] } \]
\(\sigma_D,\ f_0,\ B,\ T,\ \rho,\ \mathrm{SNR}\)
as in (11.13)

Since \(f_0^3(B^3 + 12B)\) \(= B_{\mathrm{abs}}^3 + 12 f_0^2 B_{\mathrm{abs}}\), the \(12B\) term dominates for any \(B < 2\sqrt{3}\): \(\sigma \propto 1/(f_0\sqrt{B_{\mathrm{abs}} T})\), not \(1/\sqrt{B^3 T}\). Worked: \(3/(2\pi^2 \cdot 1.25\times10^{20} \cdot 2\times10^{-6} \cdot 7.416)\) = \(8.2\times10^{-17}\) s²; bracket 0.002 at 30 dB, \(\rho\) = 1 → 0.41 ns.

Local (small-error) bound for unbiased estimators; flat coherence across a rectangular band; Gaussian statistics; no strain inside the kernel. Above threshold only.
Below a threshold SNR or correlation, false peaks dominate and the CRLB is meaningless.
knee ≈ 10 dB (2 µs, 60 %)errors in steps of \(1/f_0\)ceiling 346 ns (uniform ±24 lags)

The CRLB uses only the curvature of the main peak. When noise lifts a side peak at \(\pm m/f_0\) above the main one, the error jumps by 200 ns (154 µm at 5 MHz). The error distribution becomes a mixture, and its variance \(p\cdot(200\ \mathrm{ns})^2\) swamps the local term: at 10 dB 3–5 % of trials hop and \(\sigma\) is 10–18× the CRLB (m11_jitter_vs_crlb.m, Demo 11.3); at 5 dB a third hop. Global bounds (Ziv–Zakai, Barankin) reproduce the knee.

Approx. modelLocal bound (11.13) against the jitter floor \(\sqrt{p}\cdot 200\) ns of a hop fraction \(p\) (4 % and 33 %, the Monte Carlo values near 10 and 5 dB).
The knee moves up with short kernels (2.5 periods: hops persist at 40 dB), narrow bands, large search ranges and low \(\rho\). Report outlier fraction and a robust spread (1.48 MAD) next to \(\sigma\).
False peaks sit at multiples of the RF period because the RF correlation is a cosine under its envelope.
side peak \(\approx \exp(-\pi^2 B^2/(8\ln 2))\)0.53 at 60 %0.93 at 20 %

For a Gaussian pulse of −6 dB fractional bandwidth \(B\), the expected RF autocorrelation is \(\exp(-\tau^2/4\sigma_t^2)\cos(2\pi f_0\tau)\) with \(\sigma_t = \sqrt{2\ln 2}/(\pi B f_0)\). Every maximum of the cosine, near \(\tau = m/f_0\), is a candidate peak whose height is the envelope there: \(\exp(-\pi^2 B^2/(8\ln 2))\) at \(\tau = 1/f_0\). At 60 % that is 0.53 (the exact local maximum, pulled slightly inward by the envelope, is 0.54); at 20 % it is 0.93, so a little noise or decorrelation suffices to flip the choice.

Exact formula
Exact for the expected correlation of a Gaussian spectrum; a single finite kernel fluctuates around it. Envelope (\(\lvert\mathrm{IQ}\rvert\)) correlation has no side peaks, which is why it is used for the coarse search.
Normalize the correlation so the peak value becomes a quality metric.
\(\mathrm{E}[\hat{\rho}_{\max}] \approx \rho\,\mathrm{SNR}_p/(1+\mathrm{SNR}_p)\)\(\sigma \propto \sqrt{1/\hat{\rho}^2 - 1}\)threshold 0.7–0.9

NCC divides by the energies of the two windows, so its peak estimates the correlation coefficient between them: speckle similarity \(\rho\) times the noise coherence \(\mathrm{SNR}_p/(1+\mathrm{SNR}_p)\). The Walker–Trahey bracket is \(1/\lvert\gamma\rvert^2 - 1\), so the expected jitter is a monotone function of the peak: a peak of 0.99 predicts about 2.5× less jitter than 0.94 (\(\sqrt{0.0203/0.1317} \approx 0.39\)). Masking or weighting by peak NCC is therefore principled, not cosmetic.

Exact formula
A high peak does not prove the right lobe (side peaks can exceed 0.9 in narrowband data). Short kernels inflate the NCC of unrelated windows; the peak is a biased, noisy estimate of \(\rho\).
Decorrelation sets a jitter floor that no amount of SNR removes.
floor \(\propto \sqrt{1/\rho^2 - 1} \approx \sqrt{2(1-\rho)}\)SNR stops helping at \(1/\mathrm{SNR}_p \approx 1 - \rho\)

Expand the bracket: \((1 + 1/\mathrm{SNR}_p)^2/\rho^2 - 1\) \(\approx 2/\mathrm{SNR}_p + 2(1-\rho)\). The two terms are equal when \(\mathrm{SNR}_p = 1/(1-\rho)\): for \(\rho\) = 0.99 that is 20 dB, and at 30 dB the bound is already 1.35 ns against 0.41 ns at \(\rho\) = 1. Strain, lateral or elevational motion, and inter-element distance in phase-screen work all act through \(\rho\).

Exact formula
Treats decorrelation as an extra additive, independent component with the same spectrum as the signal; strain-induced decorrelation is structured (companding can undo much of it).
Kasai: \(v = \frac{c\cdot\mathrm{PRF}}{4\pi f_0}\angle R(1)\), and the Nyquist velocity is \(c\cdot\mathrm{PRF}/(4f_0)\).
0.385 m/s at 5 MHz, 5 kHz\(\lambda/4\) = 77 µm per pulsevar 1.8× CRLB at \(M\) = 8

A scatterer moving at \(v\) changes the two-way path by \(2v/\mathrm{PRF}\) per pulse, rotating the slow-time IQ phase by \(\Delta\varphi = 4\pi f_0 v/(c\cdot\mathrm{PRF})\) (11.14). \(R(1) = \sum z[m+1]\,z^*[m]\) averages that rotation over pulses and gates; invert for \(v\). \(\angle R(1)\) wraps at \(\pm\pi\), so \(\lvert v\rvert < c\cdot\mathrm{PRF}/(4f_0)\) = 1540·5000/(2·10⁷) = 0.385 m/s. At high SNR the phase sum telescopes to \((\varphi_{M-1} - \varphi_0)/(M-1)\): \(\sigma \propto 1/(M-1)\), equal to the CRLB at \(M\) = 2–3 and 14× above it in variance at \(M\) = 64 (20 dB).

Exact formula
Uses the nominal \(f_0\) (downshift biases \(v\) low); clutter must be filtered first; spectral broadening and low SNR add a \(1/M\) variance term; aliasing beyond \(v_{\mathrm{Nyq}}\).
Bigger kernels lower jitter but blur displacement and decorrelate under strain: there is an optimum kernel length.
1 % strain, 20 dB: \(T_{\mathrm{opt}} \approx\) 1.5 µs\(\varepsilon T_{\mathrm{opt}} f_0 \approx 0.075\)Monte Carlo 2.5 ns vs bound 2.1 ns

The noise term in the bound falls as \(1/T\), while strain decorrelation \(1 - \rho^2\) grows as \((\varepsilon T f_0)^2\), so \(\sigma(T)\) has a minimum where the two are comparable. At \(\varepsilon\) = 1 %, 20 dB, 5 MHz, 60 %: \(T_{\mathrm{opt}}\) = 1.5 µs (7.5 periods, \(\sigma\) = 2.09 ns = 1.6 µm) and axial kernel \(cT/2\) = 1.16 mm. The Monte Carlo in m11_kernel_strain.m finds its minimum at 1–1.5 µs (2.5 ns). At 3 % strain the optimum falls to 0.5 µs.

Approx. modelBound (11.13) with \(\rho(\varepsilon, T)\) from the window-averaged Gaussian-pulse correlation under uniform strain.
Uniform strain, no companding, bound-based. Short kernels (< ~5 periods here) enter the false-peak regime, so the practical optimum sits at or above the bound's. Resolution, not jitter, may set \(T\) in a heterogeneous target.
Cumulative summation of adjacent-element delays is a random walk: the error grows as \(\sqrt{n}\).
\(\sigma_n = \sqrt{n}\,\sigma_\Delta\)2 ns → 16 ns (29°) at element 64one hop = \(1/f_0\) step

Flax & O'Donnell integrate \(\Delta\hat{\tau}_i = \hat{\tau}_{i+1} - \hat{\tau}_i\). If the pair errors are independent, \(\mathrm{Var}(\hat{\tau}_n) = n\sigma_\Delta^2\) (11.17). With \(\sigma_\Delta\) = 2 ns, \(\sqrt{63} \times 2\) = 16 ns, 29° at 5 MHz. In m11_phase_screen_estimation.m (\(\rho\) = 0.984, the VCZ value for a uniform 64-element transmit) \(\sigma_\Delta\) = 1.7 ns and the far-end error is 12 ns against \(\sqrt{63}\,\sigma_\Delta\) = 13.8 ns. A single side-peak hop adds 200 ns to every later element.

SimulationThree seeded Gaussian random walks with \(\sigma_\Delta\) = 2 ns per step (JS).
Independence fails when noise is shared between pairs. Cures: beamsum reference (no accumulation), multi-lag least squares (helps when electronic noise dominates; little when speckle decorrelation dominates), NCC weighting, iteration. Constant and linear terms are unobservable.
Lateral decorrelation scales with the beam width: half a FWHM of lateral motion already costs \(\rho\) = 0.71.
\(\rho = \exp[-2\ln 2\,(\Delta x/w)^2]\)0.92 at \(w/4\)0.25 at \(w\)

Moving the scatterers laterally by \(\Delta x\) changes the speckle by the overlap of the PSF with its shifted copy. For a Gaussian pulse-echo amplitude PSF of FWHM \(w\), \(\rho = \exp(-\Delta x^2/4\sigma_x^2)\) \(= \exp[-2\ln 2\,(\Delta x/w)^2]\). With \(w\) = 0.6 mm, 0.15 mm of motion gives 0.92, which by itself caps the axial jitter at the \(\rho\) = 0.92 floor. Keep inter-frame motion small compared with the beam (raise the frame rate) or track laterally too.

Exact formula
Gaussian, shift-invariant PSF and pure translation. Elevational motion follows the slice-thickness profile, usually wider; rotation and shear decorrelate parts of the kernel differently.
Lateral tracking is an order of magnitude worse than axial, and on RF the axial peak should be fitted before the lateral one.
\(\sqrt{\text{curvature ratio}}\) 15 (RF), 2.6 (envelope)lateral 15.5 → 2.7 µm with axial-first

Jitter scales as \(1/\sqrt{-R''(0)}\). Axially the RF correlation has the carrier: \(-R'' = (2\pi f_{z0})^2 + 1/(2\sigma_z^2)\); laterally only the beam envelope: \(1/(2\sigma_x^2)\). For 5 MHz and 0.6 mm FWHM, \(\sqrt{\text{ratio}}\) = 14.9. In m11_block_matching_2d.m a separable fit on the integer axial row gives 15.5 µm lateral rms; fitting axial parabolas per column first and the lateral parabola through their peak heights gives 2.7 µm, a ratio of 13 to the axial 0.21 µm.

Exact formula
The coupling arises because the integer axial row sits up to half a sample off a carrier-modulated ridge. Transverse-oscillation beamforming or angle compounding add a lateral carrier and change the ratio.
Differentiating displacement into strain amplifies jitter by \(\sqrt{2}/\Delta z\).
\(\sigma_\varepsilon \approx \sqrt{2}\,\sigma_u/\Delta z\)0.3 µm, 1 mm → 0.04 %

A two-point difference of independent estimates has variance \(2\sigma_u^2\), divided by \(\Delta z^2\). With \(\sigma_u\) = 0.3 µm (about the 0.4 ns bound) and \(\Delta z\) = 1 mm, \(\sigma_\varepsilon\) = 0.042 % strain, so 1 % applied strain gives a strain SNR near 24 before decorrelation. This is the logic behind the strain filter: the usable strain range is bounded below by jitter and above by decorrelation.

Exact formula
Overlapping kernels make neighbouring estimates correlated (a least-squares strain estimator over several points does better); interpolator bias does not average out and adds a periodic strain artefact.

11.9Cheat sheet

ConceptFormulaUltrasound use
Two-way delay ↔ displacement\(\Delta z = cD/2\); 1 ns ↔ 0.77 µm; 1 sample (40 MHz) ↔ 19.3 µmconverting tracking output to µm
NCC\(\hat{\rho}(k) = \sum x\,y_k \big/ \sqrt{\sum x^2 \sum y_k^2}\)speckle tracking, ULM patch matching, quality mask
Expected NCC peak\(\rho\cdot\mathrm{SNR}_p/(1+\mathrm{SNR}_p)\)reject kernels with low peak NCC
Linear correlation via FFTzero-pad to \(\ge N_x + N_y - 1\)long kernels, many lags
Parabolic fit\(\delta = (R_- - R_+)/\bigl(2(R_- - 2R_0 + R_+)\bigr)\)default, bias ≈ 0.01 sample at \(f_0/f_s = 1/8\)
Cosine fit\(\delta = \arctan\bigl((R_+ - R_-)/(2R_0\sin\omega)\bigr)/\omega\)low-bias RF correlation interpolation
Phase-zero estimate\(\delta = -\angle C(k_0)/(2\pi f_0 T_s)\)IQ tracking; needs correct \(f_0\)
Bias symmetry\(\mathrm{bias}(0) = \mathrm{bias}(\tfrac{1}{2}) = 0\); \(\mathrm{bias}(d) = -\mathrm{bias}(1-d)\)pixel-locking texture in displacement maps
Walker–Trahey bound\(\sigma_D \ge \sqrt{\frac{3}{2f_0^3\pi^2 T(B^3+12B)}\bigl[\frac{(1+1/\mathrm{SNR}^2)^2}{\rho^2} - 1\bigr]}\)predicting jitter for ARFI, SWE, phase screens
Threshold / false peakserrors in steps of \(1/f_0\); dominate below ~10 dB (\(T\) = 2 µs, \(B\) = 60 %)outlier rejection, coarse-to-fine search
Kasai\(v = \frac{c\cdot\mathrm{PRF}}{4\pi f_0}\,\angle R(1)\)colour Doppler, ARFI, SWE
Nyquist velocity\(v_{\mathrm{Nyq}} = c\cdot\mathrm{PRF}/(4f_0)\)PRF choice; \(\pm\lambda/4\) per pulse
Phase-screen integration\(\sigma(\hat{\tau}_n) = \sqrt{n}\,\sigma_\Delta\)use beamsum reference or multi-lag LS to avoid drift
Kernel optimum under strainmin over \(T\) of \(\sigma_{\mathrm{WT}}(T, \rho(\varepsilon, T))\); \(\varepsilon T \approx 0.075/f_0\) at 20 dBelastography kernel design, strain filter
False-peak height (Gaussian pulse)\(\exp(-\pi^2 B^2/(8\ln 2))\): 0.53 at \(B\) = 60 %, 0.93 at 20 %predicting peak hopping, choosing coarse search
Phase-estimator range\(\lvert D\rvert < 1/(2f_0)\) (\(\pm\lambda/4\) displacement); bias \((k_0 - d)\,e/(1+e)\) for \(f_0\) error \(e\)ARFI/SWE tracking, Loupas vs Kasai
Lateral decorrelation\(\rho = \exp[-2\ln 2\,(\Delta x/w)^2]\), \(w\) = amplitude FWHMframe-rate and lateral-tracking design
Strain noise\(\sigma_\varepsilon \approx \sqrt{2}\,\sigma_u/\Delta z\)elastography SNR
Ziv–Zakai (uniform prior, width \(A\))\(\mathrm{MSE} \ge \frac{1}{A}\int_0^A h\,(A - h)\, P_{\min}(h)\, dh\)threshold analysis beyond the CRLB

11.10Interview and qualifying-exam questions

I1Why is the cross-correlation peak the maximum-likelihood delay estimate, and when is it not? · screen

For a known template in white Gaussian noise, the log-likelihood is \(-\sum (y - s_D)^2/2\sigma^2\); the energy of the shifted template does not depend on \(D\), so the likelihood is maximized by maximizing \(\sum y\,s_D\), the correlation. With both records noisy and the signal unknown, the ML processor is Knapp & Carter's generalized cross-correlation: weight the cross-spectrum by \(\lvert\gamma\rvert^2/[\lvert G_{xy}\rvert(1 - \lvert\gamma\rvert^2)]\). For white noise and high in-band SNR that is flat in band and small outside, so band-pass filtering plus plain correlation is close to ML. It fails for coloured noise or clutter (pre-whiten), for shape changes between records (strain, decorrelation), and below threshold, where the ML estimate is dominated by false peaks. Follow-up: what does PHAT weighting (\(1/\lvert G_{xy}\rvert\)) do? It whitens to phase only: sharp peaks and robustness to reverberation, but it amplifies low-SNR bands, so it is ML only when the coherence is constant across the band.

I2Compare parabolic, cosine and phase-zero subsample estimators: bias, variance, cost, failure modes. · onsite

All three start from the integer NCC peak. Parabolic: model-free and cheap, but the wrong shape for an RF peak; bias 0.012 sample at \(f_s/f_0\) = 8 and 0.052 at \(f_s/f_0\) = 4, toward the integer sample (pixel locking). Cosine: fits \(R \propto \cos(\omega(k - \delta))\), estimates \(\omega\) from the data, bias ≈ 0.0014 sample, residual from the envelope; same cost. Phase-zero: \(\delta = -\angle C(k_0)/(2\pi f_0 T_s)\) on the analytic signals, no interpolation bias if \(f_0\) equals the spectral centroid; biased by \((k_0 - d)\,e/(1+e)\) for an \(f_0\) error \(e\). Above threshold all three have jitter near the CRLB (Demo 11.3); bias, not jitter, separates them. Every symmetric fit is exact at \(d\) = 0 and ½. Follow-up: why is the Gaussian fit the worst on RF and the best on envelope data? Because it models a Gaussian peak: right for \(\lvert\mathrm{IQ}\rvert\) correlation, too sharp-shouldered for a cosine, so on RF it overshoots toward half-samples with 0.025 sample bias.

I3Derive the Kasai estimator and its Nyquist velocity. · screen

Slow-time IQ at a fixed gate: \(z[m] = A\,e^{j(\varphi_0 + m\Delta\varphi)} + n\), with the two-way path change per pulse \(2v/\mathrm{PRF}\), so \(\Delta\varphi = 2\pi f_0 (2v/c)/\mathrm{PRF}\) \(= 4\pi f_0 v/(c\cdot\mathrm{PRF})\). The lag-one autocorrelation \(R(1) = \sum z[m+1]\,z^*[m]\) \(\approx \lvert A\rvert^2 (M-1)\, e^{j\Delta\varphi}\), so \(\hat{v} = c\cdot\mathrm{PRF}\cdot\angle R(1)/(4\pi f_0)\). \(\angle R(1) \in (-\pi, \pi]\), so \(\lvert v\rvert < c\cdot\mathrm{PRF}/(4f_0)\): 0.385 m/s at 5 MHz, 5 kHz, i.e. \(\lambda/4\) = 77 µm per pulse. Equivalently \(\angle R(1)/(2\pi T_{\mathrm{PRF}})\) is the mean Doppler frequency of the slow-time spectrum. Follow-up: is it optimal? Not for a tone at high SNR: the sum telescopes to the first and last pulses, \(\sigma \propto 1/(M-1)\), variance 1.8× the CRLB at \(M\) = 8. Loupas adds the fast-time lag to estimate the local \(\bar{f}\) and remove downshift bias.

I4Why do false peaks appear at multiples of the RF period, and how do you suppress them? · onsite

The expected RF correlation is the inverse transform of the power spectrum: an envelope of width ~\(1/B\) times \(\cos(2\pi f_0\tau)\). Each cosine maximum near \(\tau = m/f_0\) is a local peak with height equal to the envelope there; for a Gaussian pulse, \(\exp(-\pi^2 B^2/(8\ln 2))\) at \(m = 1\): 0.53 at 60 % bandwidth, 0.93 at 20 %. Noise, decorrelation and the finite kernel can push one above the main peak, giving errors of exactly one period (200 ns, 154 µm at 5 MHz). Suppression: wider bandwidth; longer kernels (at least ~5 periods); a coarse search on envelope or \(\lvert\mathrm{IQ}\rvert\) correlation, then RF or phase refinement; small search ranges guided by neighbours or a previous frame; spatial regularization or median filtering of the displacement field; rejecting low-NCC kernels. Follow-up: why do hops persist at high SNR with a 2.5-period kernel? A short window of speckle is nearly a sinusoid, so the window one period away matches it almost as well; the finite-kernel fluctuation, not the noise, decides.

I5State the Walker–Trahey bound and give its scalings. What jitter do you expect at 5 MHz, 60 % bandwidth, 2 µs kernel, 30 dB? · onsite

\(\sigma_D \ge \sqrt{\frac{3}{2\pi^2 f_0^3 T(B^3+12B)}\bigl[\frac{(1+1/\mathrm{SNR}^2)^2}{\rho^2} - 1\bigr]}\) (11.13), the Carter coherence CRLB for a rectangular band with flat coherence \(\rho\cdot\mathrm{SNR}_p/(1+\mathrm{SNR}_p)\). Scalings: \(f_0^{-3/2}\), \(T^{-1/2}\), \(\approx B^{-1/2}\) for \(B < 1\), and a decorrelation floor \(\propto \sqrt{1/\rho^2 - 1}\). Numbers: prefactor \(8.2\times10^{-17}\) s², bracket 0.002 → 0.41 ns (0.31 µm); with \(\rho\) = 0.99, 1.35 ns. Follow-up: your simulation reaches 0.3 ns at 40 dB but the bound says 0.13 ns. Why? Interpolator bias and window-edge decorrelation (\(W - \lvert D\rvert\) overlap) set a floor; the bound assumes an unbiased estimator and an infinite stationary record.

I6How do the CRLB and the Ziv–Zakai bound relate, and what is the threshold? · research

The CRLB is local: it uses the Fisher information, the curvature of the log-likelihood at the true delay, so it only describes small errors around the right peak. The Ziv–Zakai bound converts estimation into binary detection: for a delay with a uniform prior of width \(A\), \(\mathrm{MSE} \ge \frac{1}{A}\int_0^A h\,(A - h)\,P_{\min}(h)\,dh\), where \(P_{\min}(h)\) is the minimum error probability of deciding between delays \(D\) and \(D + h\) (a valley-filling refinement tightens it). At high SNR only small \(h\) contribute and it converges to the CRLB. At low SNR \(P_{\min}(h)\) stays near ½ even for large \(h\), notably near \(h = m/f_0\) where the correlation side peaks are high, and the bound rises toward the prior variance \(A^2/12\) (a standard deviation of 346 ns for ±24 lags at 40 MHz). The SNR where the side-peak terms switch on is the threshold. Follow-up: how would you design around it? Increase \(T\cdot B\) or SNR until the side-peak error probability is negligible, shrink \(A\) with prior information (tracking, coarse-to-fine), and report outlier rates separately.

I7How do you choose the kernel length for quasi-static elastography or speckle tracking? · onsite

Three constraints. (1) Jitter: \(\sigma \propto T^{-1/2}\) through the bound. (2) Decorrelation: under strain \(\varepsilon\) the window ends drift by \(\varepsilon T/2\), and \(1 - \rho^2\) grows as \((\varepsilon T f_0)^2\); the optimum is where the noise and decorrelation terms balance, about \(\varepsilon T f_0 \approx 0.075\) at 20 dB (\(T \approx\) 1.5 µs at 1 %, 5 MHz, \(\sigma \approx\) 2.1 ns; Monte Carlo 2.5 ns). (3) Resolution and false peaks: axial resolution ≈ \(cT/2\) (1.16 mm for 1.5 µs), and kernels below ~5 periods hop. In practice: 5–15 periods with 50–80 % overlap, companding or multi-step (coarse-to-fine) tracking for strains above 2–3 %, and a final displacement-gradient estimator over several points. Follow-up: what is the strain filter? Varghese & Ophir's curve of elastographic SNR versus applied strain: bounded at low strain by jitter and at high strain by decorrelation.

I8How would you estimate a near-field phase screen from channel data, and what limits it? · research

Model a thin delay screen \(\tau_i\) at the aperture. Cross-correlate each element with its neighbour (Flax & O'Donnell) on speckle or a point target, integrate the \(\Delta\hat{\tau}\), apply the correction, re-transmit and iterate, since a better transmit focus raises the inter-element coherence. Limits: (1) \(\rho\) between neighbours is set by van Cittert–Zernike (\(1 - 1/N\) at focus for a uniform \(N\)-element transmit), so jitter has a decorrelation floor; (2) integration is a random walk, \(\sqrt{n}\,\sigma_\Delta\) (\(\sqrt{63} \times 2\) ns ≈ 16 ns at element 64 for \(\sigma_\Delta\) = 2 ns), and one false peak adds a \(1/f_0\) step to every later element; (3) mean and linear terms are unobservable from diffuse speckle; (4) the thin-screen, single-isoplanatic-patch assumption fails for distributed aberrators, where the correction is only valid near the region used. Better: correlate against a beamsum reference, multi-lag least squares with NCC weights, or estimate from the full set of pairs (distortion or reflection matrix methods). Follow-up: why does multi-lag LS sometimes not help? If pair errors are dominated by Markov-like speckle decorrelation, a lag-\(l\) pair is as noisy as \(l\) adjacent pairs and carries little new information; longer lags also hop more.

I9A displacement map shows a faint sawtooth texture with a period of one axial sample. What is it, and how do you remove it? · screen

Pixel (peak) locking: the subsample interpolator is biased toward integer lags, with zero bias at 0 and ½ sample and extremes in between, so a smoothly varying true displacement acquires a periodic error of one-sample period. At \(f_s/f_0\) = 8 the parabolic amplitude is ~0.012 sample (0.23 µm); after differentiation it becomes a strain artefact. Fixes: cosine fit or band-limited upsampling of the correlation, phase-zero estimation on IQ with the local \(\bar{f}\) (Loupas), or a higher sampling rate. Averaging does not help because the error is deterministic in \(d\). Follow-up: can you calibrate and subtract the bias curve? Only approximately: each kernel's interpolation error varies with its own local spectrum (±0.011 sample in m11_ncc_subsample_bias.m), comparable to the mean bias.

I10Your ARFI phase tracker underestimates displacement by about 15 % at 4 cm depth in liver, but not near the surface. Why? · onsite

Frequency-dependent attenuation shifts the received centre frequency down with depth (\(f_c = f_0 - 2\alpha_0\sigma_f^2 z\) for a Gaussian spectrum: 0.75 MHz at 4 cm for 0.5 dB/cm/MHz, 5 MHz, 60 %). The phase measures \(2\pi f_c D\), and dividing by the nominal \(f_0\) gives \(\hat{D} = D\cdot f_c/f_0 = 0.85\,D\). Fix: Loupas' 2D autocorrelator, which estimates the local mean frequency from the fast-time lag-one correlation, or a depth-dependent \(f_c(z)\) from the spectral centroid. Follow-up: does NCC with parabolic interpolation have the same problem? No scale error, because it measures time directly; its bias depends on \(f_s/f_c\) and is small at 40 MHz.

I11Why is lateral speckle tracking so much worse than axial, and what would you do about it? · onsite

Jitter scales as \(1/\sqrt{-R''(0)}\). Axially the RF correlation carries the \(f_0\) carrier, so its curvature is \((2\pi f_{z0})^2 + 1/(2\sigma_z^2)\); laterally there is only the beam envelope, \(1/(2\sigma_x^2)\). For 5 MHz and a 0.6 mm lateral FWHM the predicted ratio is about 15. Lateral sampling is also coarse (line spacing comparable to the beam), and on RF a separable 2D fit through the integer axial row couples axial subsample offsets into the lateral estimate (15.5 µm vs 2.7 µm with axial-first fitting in m11_block_matching_2d.m). Remedies: fit axially first or fit a 2D model with a cross term; denser lateral sampling or interpolated beams; transverse-oscillation beamforming or angle compounding to create a lateral carrier; combine axial estimates from several steering angles. Follow-up: why does the envelope help laterally but hurt axially? It removes the carrier: the axial curvature drops (axial error 3–4× worse) and the coupling disappears.

I12How would you use time-delay estimation to measure local sound speed? · research

Beamform the same region from several transmit (or receive) steering angles with an assumed \(c\). A local sound-speed error makes the echo from a given point arrive with angle-dependent delay errors, a line integral of the slowness error along each path. Measure the phase shift between complex images from neighbouring angles with a Kasai-type lag-one product (CUTE, Jaeger et al. 2015), stack these into a linear system relating phase shifts to the slowness map, and solve it with regularization. At 5 MHz, 1 ns of delay is 1.8° of phase, so the \(\pm\lambda/4\) wrap and the phase jitter fix the usable angle step and the averaging region; the Walker–Trahey bound with the inter-angle coherence predicts the noise. Follow-up: what limits it? Decorrelation between angles (larger steps lower \(\rho\)), the ill-conditioning of limited-angle tomography (poor axial resolution of \(c\)), and the \(f_0\) used to convert phase to delay.

11.11MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run. Runtimes are Octave timings; MATLAB is usually faster.

Download the scripts: m11_ncc_subsample_bias.m · m11_jitter_vs_crlb.m · m11_kernel_strain.m · m11_kasai_doppler.m · m11_phase_screen_estimation.m · m11_block_matching_2d.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m11_ncc_subsample_bias.mBias curves of parabolic, cosine, Gaussian and upsampled-sinc fits on noise-free shifted speckle · change fs, f0, B6 · runtime ≈ 15 s 
%% m11_ncc_subsample_bias.m
% Purpose : bias of subsample delay estimators on band-limited RF speckle.
%           Two records share the same speckle; the second is delayed by a
%           fractional shift d (0..1 sample) with an exact FFT phase shift.
%           NCC is computed at integer lags, then refined with parabolic,
%           cosine, Gaussian and upsampled-sinc fits.
% Change  : f0, fs (fs/f0 sets the parabolic bias, roughly as (f0/fs)^2),
%           B6 (-6 dB fractional bandwidth), kernel length W.
% Expect  : noise-free, so the curves are bias only. At f0 = 5 MHz,
%           fs = 40 MHz, B6 = 0.6: max |bias| about 0.012 smp (parabolic),
%           0.025 (Gaussian on RF), 0.002 (cosine), < 0.001 (upsample x8).
%           Every curve is zero at d = 0, 0.5 and 1 (within the Monte Carlo
%           standard error) and odd about d = 0.5. The noise-free 'std over
%           speckle' (0.011 smp parabolic, 0.007 cosine, 0.0008 upsample)
%           shows that each kernel's interpolation error depends on its own
%           local spectrum, so the mean bias curve cannot simply be
%           subtracted as a correction.
rng(1);

%% Parameters
c   = 1540;          % speed of sound (m/s)
f0  = 5e6;           % centre frequency (Hz)
fs  = 40e6;          % sampling rate (Hz)
B6  = 0.6;           % -6 dB fractional bandwidth of the amplitude spectrum
N   = 512;           % record length (samples)
n0  = 200;           % kernel start (sample index, 1-based)
W   = 80;            % kernel length (samples) = 2 us
K   = 24;            % search range +-K lags
L   = 12;            % half-length of the Hann-tapered sinc interpolator (lags)
U   = 8;             % upsampling factor for the sinc interpolator
nReal = 400;         % speckle realizations averaged per shift
dGrid = linspace(0, 1, 21);   % true fractional shifts (samples)

%% Band-limited speckle generator (Gaussian amplitude spectrum)
f   = (0:N-1)' * fs / N;  f(f >= fs/2) = f(f >= fs/2) - fs;   % signed bin frequencies (Hz)
sf  = B6 * f0 / 2 / sqrt(2*log(2));                       % spectral std (Hz)
H   = exp(-(abs(f) - f0).^2 / (2*sf^2));                  % even amplitude response
H(N/2+1) = 0;                                             % no energy at Nyquist
snc = @(x) (sin(pi*x) + (x == 0)) ./ (pi*x + (x == 0));   % sinc without toolbox

%% Estimators applied to three NCC samples Rm, R0, Rp around the integer peak
fitPar = @(Rm, R0, Rp) 0.5 * (Rm - Rp) ./ (Rm - 2*R0 + Rp);
fitCos = @(Rm, R0, Rp) atan((Rp - Rm) ./ (2*R0 .* sin(acos((Rm + Rp) ./ (2*R0))))) ./ acos((Rm + Rp) ./ (2*R0));
fitGau = @(Rm, R0, Rp) 0.5 * (log(Rm) - log(Rp)) ./ (log(Rm) - 2*log(R0) + log(Rp));

lags = -K:K;
est  = zeros(numel(dGrid), nReal, 4);     % parabolic, cosine, Gaussian, upsampled sinc
idxW = (n0:n0+W-1)';  Yidx = bsxfun(@plus, idxW, lags);
for r = 1:nReal
    S  = fft(randn(N, 1)) .* H;           % speckle spectrum
    x  = real(ifft(S));
    for id = 1:numel(dGrid)
        d = dGrid(id);
        y = real(ifft(S .* exp(-1j*2*pi*f*d/fs)));   % exact fractional delay
        % NCC: fixed kernel of x, sliding window of y
        xw = x(idxW);  Yw = y(Yidx);                       % W x (2K+1) sliding windows
        R  = (xw' * Yw) ./ sqrt(sum(xw.^2) * sum(Yw.^2, 1));
        [~, im] = max(R(2:end-1));  im = im + 1;  k0 = lags(im);
        Rm = R(im-1); R0 = R(im); Rp = R(im+1);
        est(id, r, 1) = k0 + fitPar(Rm, R0, Rp);
        est(id, r, 2) = k0 + fitCos(Rm, R0, Rp);
        if min([Rm R0 Rp]) > 0, est(id, r, 3) = k0 + fitGau(Rm, R0, Rp); else, est(id, r, 3) = NaN; end
        % band-limited (sinc) interpolation of R on a 1/U grid, parabolic fit on the fine grid
        tf = k0 + (-U:U) / U;                              % fine lag grid (samples)
        u  = bsxfun(@minus, tf, lags');                    % lag offsets (numel(lags) x numel(tf))
        Rf = R * (snc(u) .* (0.5 + 0.5*cos(pi*u/(L+1))) .* (abs(u) <= L));
        [~, j] = max(Rf(2:end-1));  j = j + 1;
        est(id, r, 4) = tf(j) + fitPar(Rf(j-1), Rf(j), Rf(j+1)) / U;
    end
end

%% Bias = mean error over realizations; jitter from finite-window speckle
err  = bsxfun(@minus, est, dGrid');                       % (shift, realization, method)
bias = zeros(numel(dGrid), 4); jit = zeros(1, 4);   % jit: spread over speckle (finite kernel)
for m = 1:4
    e = err(:, :, m);
    for id = 1:numel(dGrid), v = e(id, :); bias(id, m) = mean(v(~isnan(v))); end
    v = e(:); v = v(~isnan(v)); jit(m) = std(v);
end
names = {'parabolic', 'cosine', 'Gaussian', 'upsample x8'};
fprintf('f0 = %.1f MHz, fs = %.1f MHz (fs/f0 = %.1f), B6 = %.0f %%, W = %d smp\n', f0/1e6, fs/1e6, fs/f0, 100*B6, W);
fprintf('1 sample = %.1f ns = %.2f um axial\n', 1e9/fs, c/(2*fs)*1e6);
for m = 1:4
    [mx, i] = max(abs(bias(:, m)));
    fprintf('%-12s max|bias| = %.4f smp (%.3f ns) at d = %.2f | std over speckle = %.4f smp\n', ...
        names{m}, mx, mx*1e9/fs, dGrid(i), jit(m));
end
fprintf('bias at d = 0 / 0.5 (parabolic): %.1e / %.1e smp (Monte Carlo standard error %.1e)\n', ...
    bias(1,1), bias(11,1), jit(1)/sqrt(nReal));

%% Figure
figure;
subplot(1, 2, 1);
plot(dGrid, bias(:,1), 'o-', dGrid, bias(:,2), 's-', dGrid, bias(:,3), '^-', dGrid, bias(:,4), 'd-');
hold on; plot([0 1], [0 0], 'k:'); hold off;
xlabel('true fractional shift d (samples)'); ylabel('bias (samples)');
title('Bias = mean estimate - true'); legend(names, 'Location', 'southoutside'); grid on;
subplot(1, 2, 2);
lagAx = (-24:24); tt = lagAx / fs; st = 1 / (2*pi*sf);
Rth = exp(-tt.^2 / (4*st^2)) .* cos(2*pi*f0*tt);          % expected RF autocorrelation
tfine = linspace(-24, 24, 2000) / fs;
plot(tfine*1e9, exp(-tfine.^2/(4*st^2)) .* cos(2*pi*f0*tfine), 'Color', [0.6 0.6 0.6]); hold on;
plot(tt*1e9, Rth, 'o'); hold off;
xlabel('lag (ns)'); ylabel('expected NCC'); title('Expected RF NCC'); xlim([-400 400]); grid on;

What to look for. At 5 MHz, 40 MHz, 60 %: max |bias| 0.0123 sample (0.31 ns) parabolic at \(d\) = 0.30, 0.0017 cosine, 0.0251 Gaussian (opposite sign), 0.0003 upsample ×8, matching Demo 11.2. The bias at \(d\) = 0 and ½ is zero within the Monte Carlo standard error (\(5.5\times10^{-4}\)). The noise-free "std over speckle" column (0.011 sample for the parabola) is the kernel-to-kernel variation of the interpolation error. Set fs = 20e6 to see the parabolic bias grow about 4×.

m11_jitter_vs_crlb.mMonte Carlo σ of NCC + parabola and phase estimators vs SNR, against the Carter CRLB and the Walker–Trahey form · change W, B6, rho · runtime ≈ 2 s 
%% m11_jitter_vs_crlb.m
% Purpose : Monte Carlo jitter of two delay estimators versus SNR, compared
%           with the Cramer-Rao bound (Carter coherence integral for this
%           spectrum) and the Walker-Trahey closed form; outlier fraction
%           shows the threshold effect (jumps to a side peak 1/f0 away).
% Change  : W (kernel length), B6 (bandwidth), rho (speckle correlation),
%           nTrial. A shorter kernel or narrower band moves the knee to
%           higher SNR; rho < 1 adds a high-SNR floor.
% Expect  : f0 = 5 MHz, B6 = 0.6, T = 2 us, rho = 1. From 15 to 30 dB both
%           estimators sit within about 1.3x of the Carter CRLB (0.41 ns
%           at 30 dB). At 10 dB a few % of trials hop to a side peak and the
%           std is about 10x the bound; at 5 dB about a third hop. Above
%           35 dB a floor near 0.3 ns appears (parabolic bias plus window-
%           edge decorrelation). For this Gaussian spectrum the Carter
%           integral and the WT closed form evaluated at B = B6 agree to
%           within 5 %.
rng(1);

%% Parameters
f0 = 5e6;  fs = 40e6;  B6 = 0.6;      % centre frequency, sampling rate (Hz), -6 dB fractional BW
N  = 512;  n0 = 200;   W = 80;        % record, kernel start (1-based), kernel length (samples)
K  = 24;                               % search range +-K lags
rho = 1;                               % speckle correlation between the records
snrdB  = -10:5:40;                     % SNR inside the -6 dB band (dB)
nTrial = 300;                          % trials per SNR
T = W / fs;                            % kernel duration (s)

%% Spectrum (one-sided analytic construction gives RF and IQ together)
f  = (0:N-1)' * fs / N;  pos = (f > 0) & (f < fs/2);   % positive-frequency bins
sf = B6 * f0 / 2 / sqrt(2*log(2));                       % spectral std of the amplitude response (Hz)
H  = exp(-(f - f0).^2 / (2*sf^2)) .* pos;                % amplitude response, positive bins only
Hn = H / sqrt(sum(H.^2) / N^2 * 2);                      % scale so the real part has unit power
nb = sum(pos & abs(f - f0) <= B6*f0/2);                  % bins inside the -6 dB band
fitPar = @(Rm, R0, Rp) 0.5 * (Rm - Rp) ./ (Rm - 2*R0 + Rp);
idxW = (n0:n0+W-1)';  lags = -K:K;
Yidx = bsxfun(@plus, idxW, lags);                        % W x (2K+1) window indices

sdPar = zeros(size(snrdB)); sdPha = sdPar; outFrac = sdPar; madPar = sdPar;
for is = 1:numel(snrdB)
    snr = 10^(snrdB(is)/10);
    sn  = sqrt(N^2 / (2 * nb * snr));                    % noise std per positive bin: in-band real-part noise power = 1/snr
    ePar = zeros(nTrial, 1); ePha = ePar;
    for t = 1:nTrial
        d  = rand;                                       % true shift, uniform in [0,1) samples
        S  = (randn(N,1) + 1j*randn(N,1)) / sqrt(2) .* Hn;
        U  = (randn(N,1) + 1j*randn(N,1)) / sqrt(2) .* Hn;
        Sy = (rho*S + sqrt(1 - rho^2)*U) .* exp(-1j*2*pi*f*d/fs);
        n1 = (randn(N,1) + 1j*randn(N,1)) / sqrt(2) * sn .* pos;
        n2 = (randn(N,1) + 1j*randn(N,1)) / sqrt(2) * sn .* pos;
        xa = ifft(2 * (S  + n1));  ya = ifft(2 * (Sy + n2));   % analytic signals; RF = real part
        x  = real(xa);  y = real(ya);
        xw = x(idxW);  Yw = y(Yidx);
        R  = (xw' * Yw) ./ sqrt(sum(xw.^2) * sum(Yw.^2, 1));   % NCC at all lags
        [~, im] = max(R(2:end-1));  im = im + 1;  k0 = lags(im);
        ePar(t) = k0 + fitPar(R(im-1), R(im), R(im+1)) - d;
        Cc = sum(ya(idxW + k0) .* conj(xa(idxW)));              % complex correlation at k0
        ePha(t) = k0 - angle(Cc) * fs / (2*pi*f0) - d;
    end
    ePar = ePar / fs;  ePha = ePha / fs;                  % seconds
    sdPar(is) = std(ePar);  sdPha(is) = std(ePha);
    madPar(is) = 1.4826 * median(abs(ePar - median(ePar)));
    outFrac(is) = mean(abs(ePar) > 0.5 / f0);
end

%% Bounds
% Carter CRLB: var >= [2T int (2 pi f)^2 g^2/(1-g^2) df]^-1, g = rho*S(f)/(S(f)+N0)
fg = linspace(1, fs/2, 4000);  dfg = fg(2) - fg(1);
Sg = exp(-(fg - f0).^2 / sf^2);  Sg = Sg / (sum(Sg) * dfg);     % unit-power signal PSD (one-sided)
crlb = zeros(size(snrdB)); wtb = crlb;
for is = 1:numel(snrdB)
    snr = 10^(snrdB(is)/10);  N0 = 1 / (snr * B6 * f0);         % white noise PSD: in-band noise power = 1/snr
    g = rho * Sg ./ (Sg + N0);
    crlb(is) = 1 / sqrt(2 * T * sum((2*pi*fg).^2 .* g.^2 ./ (1 - g.^2 + eps)) * dfg);
    wtb(is) = sqrt(3 / (2 * f0^3 * pi^2 * T * (B6^3 + 12*B6)) * ((1 + 1/snr)^2 / rho^2 - 1));
end

%% Print
fprintf('f0 = %.1f MHz, B6 = %.0f %%, T = %.2f us (%d smp), rho = %.3f, %d trials/SNR\n', f0/1e6, 100*B6, T*1e6, W, rho, nTrial);
fprintf(' SNR(dB)  sd_par(ns)  sd_phase(ns)  CRLB(ns)  WT(ns)  sd_par/CRLB  outliers\n');
for is = 1:numel(snrdB)
    fprintf('  %4d   %9.3f   %10.3f   %8.3f  %6.3f   %8.2f    %5.1f %%\n', snrdB(is), sdPar(is)*1e9, sdPha(is)*1e9, ...
        crlb(is)*1e9, wtb(is)*1e9, sdPar(is)/crlb(is), 100*outFrac(is));
end
fprintf('uniform over +-K lags would give sd = %.1f ns\n', 2*K/fs/sqrt(12)*1e9);

%% Figure
figure;
subplot(1, 2, 1);
semilogy(snrdB, sdPar*1e9, 'o-', snrdB, sdPha*1e9, 's--', snrdB, crlb*1e9, 'k-', snrdB, wtb*1e9, 'k:', snrdB, madPar*1e9, 'x');
xlabel('in-band SNR (dB)'); ylabel('std of delay error (ns)'); grid on;
legend('NCC+parab.', 'phase', 'CRLB', 'WT (B = B6)', '1.48 MAD', 'Location', 'northeast');
title('Jitter vs SNR');
subplot(1, 2, 2);
bar(snrdB, 100*outFrac); set(gca, 'XTick', -10:10:40); xlabel('SNR (dB)'); ylabel('outliers (%)'); grid on;
title('Outliers: |error| > 1/(2 f0)');

What to look for. \(\sigma_{\mathrm{par}}/\mathrm{CRLB}\) = 1.21, 1.04, 1.20, 1.27 at 15–30 dB (CRLB 0.405 ns at 30 dB); at 10 dB 2.7 % outliers and 9.6× the bound; at 5 dB 36 % outliers; below 0 dB \(\sigma\) approaches the uniform-search ceiling of 346 ns. Above 35 dB a floor near 0.3 ns appears. The Carter integral and the WT closed form agree within 5 % for this spectrum at \(\rho\) = 1.

m11_kernel_strain.mMonte Carlo jitter vs kernel length under 1 % axial strain, against the WT bound with ρ(ε, T) · change epsStrain, snrdB · runtime ≈ 22 s 
%% m11_kernel_strain.m
% Purpose : Monte Carlo check of the kernel-length tradeoff under axial
%           strain (Demo 11.5). Point scatterers are imaged with a Gaussian
%           RF pulse, then compressed by a uniform strain eps about the
%           kernel centre and shifted by a random subsample delay d. NCC with
%           a cosine subsample fit estimates the delay at the kernel centre
%           for a range of kernel lengths T. The rms error is compared with
%           the Walker-Trahey bound evaluated with rho(eps, T), the
%           window-averaged RF correlation under strain.
% Change  : epsStrain (0 gives a monotone T^-1/2 curve), snrdB, f0, B6.
% Expect  : f0 = 5 MHz, B6 = 0.6, 20 dB in-band SNR, eps = 1 %: the bound
%           has its minimum at T = 1.5 us (2.09 ns = 1.6 um). The Monte Carlo
%           rms error has its minimum at T = 1-1.5 us (about 2.5 ns, 1.2x the
%           bound) and rises to 3.9 ns at 5 us as strain decorrelates the
%           kernel. Below 1 us side-peak hops (0.5-7 % of trials) make the
%           rms error 5-15x the bound while the robust spread (1.48 MAD)
%           still follows it: the threshold effect. With eps = 0 the error
%           falls as T^-1/2 (0.80 ns at 5 us, bound 0.81 ns).
rng(1);

%% Parameters
c  = 1540;  f0 = 5e6;  fs = 40e6;  B6 = 0.6;     % (m/s), (Hz), (Hz), -6 dB fractional BW
epsStrain = 0.01;                                 % axial strain (compression), 1 %
snrdB = 20;                                       % SNR inside the band (signal/noise power)
Tlist = [0.25 0.5 0.75 1 1.5 2 3 4 5] * 1e-6;     % kernel lengths (s)
nTrial = 200;                                     % speckle realizations
K = 12;                                           % search +-K lags
N = 640;  tc = 320;                               % record length and kernel centre (samples, 1-based)
dens = 40;                                        % scatterers per microsecond of record

%% Pulse and helpers
sigT = sqrt(2*log(2)) / (pi*B6*f0);               % Gaussian envelope std (s) for -6 dB BW = B6*f0
pulse = @(t) exp(-t.^2 / (2*sigT^2)) .* cos(2*pi*f0*t);
t = ((1:N)' - tc) / fs;                           % time axis, 0 at the kernel centre (s)
Ns = round(dens * N / fs * 1e6);                  % scatterers per record
fitCos = @(Rm, R0, Rp) atan((Rp - Rm) ./ (2*R0 .* sin(acos((Rm + Rp) ./ (2*R0))))) ./ acos((Rm + Rp) ./ (2*R0));
fitPar = @(Rm, R0, Rp) 0.5 * (Rm - Rp) ./ (Rm - 2*R0 + Rp);
lags = -K:K;
% band-limited noise: white noise filtered by the same pulse, so the coherence is flat across the band
hn = pulse((-60:60)' / fs);
err = zeros(nTrial, numel(Tlist));

for it = 1:nTrial
    tj = (rand(1, Ns) - 0.5) * N / fs;            % scatterer times (s), centred on the kernel
    aj = randn(1, Ns);                            % scattering amplitudes
    d  = rand / fs;                               % true delay at the kernel centre (s)
    x  = pulse(bsxfun(@minus, t, tj)) * aj';
    y  = pulse(bsxfun(@minus, t, tj*(1 - epsStrain) + d)) * aj';   % compressed about tc, then delayed
    sc = sqrt(mean(x.^2));
    n1 = conv(randn(N + 120, 1), hn, 'valid');  n1 = n1(1:N);
    n2 = conv(randn(N + 120, 1), hn, 'valid');  n2 = n2(1:N);
    gn = sc * 10^(-snrdB/20) / sqrt(mean([n1; n2].^2));
    x = x + gn*n1;  y = y + gn*n2;
    for iT = 1:numel(Tlist)
        W = round(Tlist(iT) * fs);
        idxW = (tc - floor(W/2) : tc - floor(W/2) + W - 1)';
        Yw = y(bsxfun(@plus, idxW, lags));  xw = x(idxW);
        R = (xw' * Yw) ./ sqrt(sum(xw.^2) * sum(Yw.^2, 1));
        [~, im] = max(R(2:end-1));  im = im + 1;
        dl = fitCos(R(im-1), R(im), R(im+1));
        if ~isfinite(dl) || abs(dl) > 1, dl = fitPar(R(im-1), R(im), R(im+1)); end
        err(it, iT) = (lags(im) + dl) / fs - d;   % kernel centre delay is d (strain is about tc)
    end
end

%% Bound: Walker-Trahey with rho(eps, T) = mean over the window of the RF ACF at lag eps*t
acf = @(tau) exp(-tau.^2 / (4*sigT^2)) .* cos(2*pi*f0*tau);
snr = 10^(snrdB/10);
Tf = linspace(0.1e-6, 6e-6, 300);  wt = zeros(size(Tf));
for i = 1:numel(Tf)
    tt = linspace(-Tf(i)/2, Tf(i)/2, 201);  rho = max(mean(acf(epsStrain * tt)), 1e-3);
    wt(i) = sqrt(3 / (2*f0^3*pi^2*Tf(i)*(B6^3 + 12*B6)) * ((1 + 1/snr)^2 / rho^2 - 1));
end
[wmin, iw] = min(wt);
rmsE = sqrt(mean(err.^2, 1));
hop  = mean(abs(err) > 0.5/f0, 1);
madE = 1.4826 * median(abs(err - median(err, 1)), 1);

%% Print
fprintf('f0 = %.1f MHz, B6 = %.0f %%, strain = %.1f %%, in-band SNR = %d dB, %d trials\n', f0/1e6, 100*B6, 100*epsStrain, snrdB, nTrial);
fprintf('bound minimum: T = %.2f us, sigma = %.2f ns (%.2f um)\n', Tf(iw)*1e6, wmin*1e9, wmin*c/2*1e6);
fprintf('  T(us)  periods  rms err(ns)  1.48 MAD(ns)  WT bound(ns)  hops\n');
for iT = 1:numel(Tlist)
    tt = linspace(-Tlist(iT)/2, Tlist(iT)/2, 201);  rho = mean(acf(epsStrain * tt));
    b = sqrt(3 / (2*f0^3*pi^2*Tlist(iT)*(B6^3 + 12*B6)) * ((1 + 1/snr)^2 / rho^2 - 1));
    fprintf('  %4.2f   %5.2f    %9.2f    %9.2f    %9.2f    %5.1f %%\n', Tlist(iT)*1e6, Tlist(iT)*f0, rmsE(iT)*1e9, madE(iT)*1e9, b*1e9, 100*hop(iT));
end

%% Figure
figure;
loglog(Tlist*1e6, rmsE*1e9, 'o-', Tlist*1e6, madE*1e9, 's--', Tf*1e6, wt*1e9, 'k-');
hold on; loglog([2 2]/f0*1e6, [0.1 1000], 'k:'); hold off;
xlabel('kernel length T (us)'); ylabel('delay error (ns)'); grid on; axis([0.1 6 0.3 300]);
legend('rms (Monte Carlo)', '1.48 MAD', 'WT bound with rho(eps,T)', 'T = 2/f0', 'Location', 'northeast');
title(sprintf('Kernel length vs strain: eps = %.1f %%, SNR = %d dB', 100*epsStrain, snrdB));

What to look for. The bound minimum is at \(T\) = 1.48 µs (2.09 ns), as in Demo 11.5. The simulated rms error is lowest at 1–1.5 µs (2.5 ns) and climbs to 3.9 ns at 5 µs. At 0.25–0.75 µs the rms error is 15–54 ns because 0.5–7 % of trials hop, while the robust 1.48 MAD (2.4–3.7 ns) still tracks the bound. Set epsStrain = 0 and the error falls monotonically as \(T^{-1/2}\) (0.80 ns at 5 µs).

m11_kasai_doppler.mKasai estimate vs true velocity with aliasing; velocity σ vs ensemble length against a high-SNR formula and the tone CRLB · change PRF, f0, M · runtime ≈ 1.5 s 
%% m11_kasai_doppler.m
% Purpose : lag-one autocorrelation (Kasai) velocity estimator on a slow-time
%           IQ ensemble: (1) estimate vs true velocity with aliasing beyond
%           v_Nyq = c*PRF/(4 f0); (2) velocity std vs ensemble length M,
%           compared with a high-SNR formula and with the CRLB for a tone.
% Change  : PRF, f0, SNR, M. Doubling PRF doubles v_Nyq and, at fixed phase
%           noise, doubles the velocity std.
% Expect  : f0 = 5 MHz, PRF = 5 kHz: v_Nyq = 0.385 m/s, lambda/4 = 77 um per
%           pulse. Part 2 (one gate, constant amplitude, SNR 20 dB): std
%           follows var(dphi) = 1/(SNR (M-1)^2) + 1/(2 SNR^2 (M-1)), so it falls
%           as 1/(M-1) when signal x noise terms dominate; it equals the CRLB
%           6/(SNR M (M^2-1)) at M = 2 and 3 and is about 1.8x the CRLB in
%           variance at M = 8: at high SNR the phase sum telescopes, so the
%           estimate effectively uses only the first and last pulses.
rng(1);

%% Parameters
c   = 1540;     % speed of sound (m/s)
f0  = 5e6;      % transmit centre frequency (Hz)
PRF = 5e3;      % pulse repetition frequency (Hz)
G   = 8;        % axial gates averaged in R(1)
M   = 8;        % ensemble length (pulses)
snrdB = 20;     % per-sample SNR (signal power / noise power)
vNyq = c * PRF / (4 * f0);
fprintf('v_Nyq = c*PRF/(4 f0) = %.3f m/s; lambda/4 = %.1f um; ensemble duration = %.2f ms\n', ...
    vNyq, c/f0/4*1e6, M/PRF*1e3);

% Kasai estimate from an ensemble z (M x G): R(1) summed over pulses and gates
kasai = @(z) c * PRF / (4*pi*f0) * angle(sum(sum(z(2:end, :) .* conj(z(1:end-1, :)))));
% Simulate a G-gate ensemble: random complex speckle amplitude per gate, constant over slow time
ens = @(v, Mm, Gg, snr) bsxfun(@times, (randn(1,Gg) + 1j*randn(1,Gg))/sqrt(2), ...
      exp(1j * (0:Mm-1)' * 4*pi*f0*v/(c*PRF))) + (randn(Mm,Gg) + 1j*randn(Mm,Gg)) / sqrt(2*snr);

%% Part 1: estimate vs true velocity (aliasing)
snr = 10^(snrdB/10);
vTrue = linspace(-1.5, 1.5, 241);
vHat = zeros(size(vTrue));
for i = 1:numel(vTrue), vHat(i) = kasai(ens(vTrue(i), M, G, snr)); end
vWrap = vTrue - 2*vNyq * round(vTrue / (2*vNyq));
fprintf('Part 1: max |v_hat - wrapped v| inside the band = %.3f m/s\n', max(abs(vHat(abs(vTrue) < 0.9*vNyq) - vWrap(abs(vTrue) < 0.9*vNyq))));
v1 = 0.5;  fprintf('v = %.2f m/s (1.3 v_Nyq) reads as %.3f m/s (expected wrap %.3f)\n', v1, kasai(ens(v1, M, G, snr)), v1 - 2*vNyq);

%% Part 2: std vs ensemble length, one gate, constant amplitude (a tone in noise)
Ms = [2 3 4 6 8 12 16 24 32 48 64];  nT = 4000;  v0 = 0.1;
sdMC = zeros(size(Ms));  sdTh = sdMC;  sdCR = sdMC;
k = c * PRF / (4*pi*f0);                         % m/s per radian
for im = 1:numel(Ms)
    Mm = Ms(im);  m = (0:Mm-1)';
    ph = 2*pi*rand(1, nT);                       % random start phase per trial
    z  = exp(1j * bsxfun(@plus, m * 4*pi*f0*v0/(c*PRF), ph)) + (randn(Mm,nT) + 1j*randn(Mm,nT)) / sqrt(2*snr);
    est = k * angle(sum(z(2:end,:) .* conj(z(1:end-1,:)), 1));
    sdMC(im) = std(est);
    sdTh(im) = k * sqrt(1/(snr*(Mm-1)^2) + 1/(2*snr^2*(Mm-1)));
    sdCR(im) = k * sqrt(6 / (snr * Mm * (Mm^2 - 1)));
end
fprintf('Part 2 (SNR %d dB, one gate, %d trials):\n   M   T_ens(ms)  sd_MC(mm/s)  sd_formula  sd_CRLB  var ratio MC/CRLB\n', snrdB, nT);
for im = 1:numel(Ms)
    fprintf('  %2d   %6.2f     %8.3f    %8.3f   %7.3f    %5.2f\n', Ms(im), Ms(im)/PRF*1e3, 1e3*sdMC(im), 1e3*sdTh(im), 1e3*sdCR(im), (sdMC(im)/sdCR(im))^2);
end

%% Figure
figure;
subplot(1, 2, 1);
vW = vWrap;  vW([false, abs(diff(vWrap)) > vNyq]) = NaN;   % break the sawtooth at the wraps
plot(vTrue, vHat, '.', vTrue, vW, '-'); hold on;
plot([vNyq vNyq], [-1.5 1.5], 'k:', [-vNyq -vNyq], [-1.5 1.5], 'k:'); hold off;
xlabel('true axial velocity (m/s)'); ylabel('Kasai estimate (m/s)'); grid on;
title(sprintf('Aliasing at v_{Nyq} = %.3f m/s', vNyq)); axis([-1.5 1.5 -0.5 0.5]);
subplot(1, 2, 2);
loglog(Ms, 1e3*sdMC, 'o', Ms, 1e3*sdTh, '-', Ms, 1e3*sdCR, 'k--');
xlabel('ensemble length M (pulses)'); ylabel('std of v (mm/s)'); grid on;
legend('Monte Carlo', 'high-SNR formula', 'CRLB (tone)', 'Location', 'southwest');
title(sprintf('Std vs M, SNR = %d dB', snrdB));

What to look for. \(v_{\mathrm{Nyq}}\) = 0.385 m/s; 0.5 m/s reads as −0.270 m/s. In part 2 the Monte Carlo \(\sigma\) follows \(\mathrm{var}(\Delta\varphi) = 1/(\mathrm{SNR}(M-1)^2)\) \(+\, 1/(2\,\mathrm{SNR}^2 (M-1))\); the variance ratio to the CRLB is 1.01 and 1.03 at \(M\) = 2, 3, then 1.83 at \(M\) = 8 and 14.3 at \(M\) = 64: the telescoping sum discards the middle pulses.

m11_phase_screen_estimation.mAdjacent-element NCC delays on a 64-element aperture, cumulative sum vs weighted multi-lag least squares · change rho, snrdB, nLagLS · runtime ≈ 11 s 
%% m11_phase_screen_estimation.m
% Purpose : near-field phase-screen (time-delay) estimation from diffuse
%           speckle on a 64-element aperture. Each element sees partially
%           correlated speckle (adjacent correlation rho) delayed by a smooth
%           random screen tau_i. Adjacent-element NCC delays are integrated
%           with a cumulative sum (Flax & O'Donnell style) and, for
%           comparison, with a least-squares fit to lags 1..3.
% Change  : rho, snrdB, screen rms and correlation length, nLagLS.
% Expect  : rho = 0.984, 20 dB: adjacent-pair jitter sigma_D = 1.7 ns,
%           close to the decorrelation-only Walker-Trahey floor (1.6 ns for
%           1/rho^2 - 1 = 0.033, T = 2 us). The cumulative-sum error grows as
%           sqrt(n)*sigma_D (a random walk): 12 ns rms (22 deg at 5 MHz) at
%           element 64 against sqrt(63)*sigma_D = 13.8 ns. Weighted least
%           squares over lags 1..3 improves this by only about 5 %: with the
%           Markov (rho^|lag|) inter-element correlation used here, a lag-l
%           pair is as noisy as l adjacent pairs and adds little independent
%           information. With rho = 0.95, about 1 % of lag-3 pairs hop by
%           1/f0 and the LS result (about 49 ns) is worse than the plain
%           cumulative sum (23 ns). When electronic noise dominates instead
%           (rho = 0.9999, snrdB = 10) LS over lags 1..3 wins clearly
%           (about 3 ns vs 7 ns). Noise is white over the full 0-20 MHz band
%           (no band-pass prefilter), which inflates sigma_D at low SNR.
rng(1);

%% Parameters
f0 = 5e6;  fs = 40e6;  B6 = 0.6;   % centre frequency, sampling (Hz), -6 dB fractional BW
Nel = 64;                          % elements (pitch 0.3 mm)
N  = 1024;  n0 = 400;  W = 80;     % record, kernel start (1-based), kernel 2 us
K  = 24;                           % search +-K lags (+-600 ns)
rho = 0.984;                       % adjacent-element speckle correlation (VCZ at focus, uniform 64-el. transmit: 1 - 1/64)
snrdB = 20;                        % channel SNR (dB, real-part power ratio over the record)
rmsNs = 30;  Lcorr = 6;            % screen rms (ns) and correlation length (elements)
nReal = 100;                       % realizations for the error statistics
nLagLS = 3;                        % lags used by the least-squares integrator

f  = (0:N-1)' * fs / N;  f(f >= fs/2) = f(f >= fs/2) - fs;          % signed frequencies (Hz)
sf = B6 * f0 / 2 / sqrt(2*log(2));
H  = exp(-(abs(f) - f0).^2 / (2*sf^2));  H(N/2+1) = 0;
fitCos = @(Rm, R0, Rp) atan((Rp - Rm) ./ (2*R0 .* sin(acos((Rm + Rp) ./ (2*R0))))) ./ acos((Rm + Rp) ./ (2*R0));
idxW = (n0:n0+W-1)';  lags = -K:K;  Yidx = bsxfun(@plus, idxW, lags);
el = (0:Nel-1)';
% Gaussian-correlated screen: filter white noise along the aperture
gk = exp(-(-3*Lcorr:3*Lcorr).^2 / (2*(Lcorr/2)^2));  gk = gk / sqrt(sum(gk.^2));
% Least-squares design matrix: tau(i+l) - tau(i) = D(i,l), tau(1) = 0
rows = [];  for l = 1:nLagLS, for i = 1:Nel-l, rows = [rows; i, i+l, l]; end, end %#ok<AGROW>
A = zeros(size(rows,1), Nel);
for r = 1:size(rows,1), A(r, rows(r,1)) = -1; A(r, rows(r,2)) = 1; end
A = A(:, 2:end);                                                     % fix tau(1) = 0

errCS = zeros(Nel, nReal);  errLS = errCS;  dAdj = zeros(Nel-1, nReal);  trueAdj = dAdj;
for r = 1:nReal
    tau = conv(randn(Nel + 6*Lcorr, 1), gk', 'valid');  tau = tau(1:Nel);
    tau = (tau - mean(tau)) / std(tau) * rmsNs * 1e-9;               % screen (s)
    % partially correlated speckle across elements: AR(1) in element index
    S = zeros(N, Nel);  S(:,1) = randn(N,1);
    for i = 2:Nel, S(:,i) = rho*S(:,i-1) + sqrt(1 - rho^2)*randn(N,1); end
    X = real(ifft(bsxfun(@times, fft(S), H) .* exp(-1j*2*pi*f*tau')));  % delayed band-limited channels
    X = X / sqrt(mean(X(:).^2));
    X = X + randn(N, Nel) * 10^(-snrdB/20);
    D = zeros(Nel, nLagLS);  Rpk = zeros(Nel, nLagLS);
    for l = 1:nLagLS
        for i = 1:Nel-l
            xw = X(idxW, i);  Yw = X(:, i+l);  Yw = Yw(Yidx);
            R = (xw' * Yw) ./ sqrt(sum(xw.^2) * sum(Yw.^2, 1));
            [~, im] = max(R(2:end-1));  im = im + 1;
            dl = fitCos(R(im-1), R(im), R(im+1));  if ~isfinite(dl), dl = 0; end
            D(i, l) = (lags(im) + dl) / fs;                         % estimated tau(i+l) - tau(i)
            Rpk(i, l) = R(im);                                      % peak NCC (quality)
        end
    end
    tauCS = [0; cumsum(D(1:Nel-1, 1))];                              % cumulative sum of adjacent delays
    % weighted least squares: weight = Fisher-information proxy rho^2/(1-rho^2) from the peak NCC,
    % then one pass of outlier rejection (residual > 1/(4 f0)) to drop side-peak hops
    b = zeros(size(rows,1), 1);  wq = b;
    for q = 1:size(rows,1), b(q) = D(rows(q,1), rows(q,3)); rq = min(Rpk(rows(q,1), rows(q,3)), 0.999); wq(q) = rq^2/(1 - rq^2); end
    sw = sqrt(wq);  t1 = (bsxfun(@times, A, sw)) \ (b .* sw);
    keep = abs(A*t1 - b) < 0.25/f0;
    if rank(A(keep,:)) < Nel-1, keep = true(size(b)); end          % never drop an element entirely
    t1 = (bsxfun(@times, A(keep,:), sw(keep))) \ (b(keep) .* sw(keep));
    tauLS = [0; t1];
    errCS(:, r) = tauCS - (tau - tau(1));
    errLS(:, r) = tauLS - (tau - tau(1));
    dAdj(:, r) = D(1:Nel-1, 1);  trueAdj(:, r) = diff(tau);
end

%% Statistics
eAdj = dAdj(:) - trueAdj(:);
hop = abs(eAdj) > 0.5/f0;
sigD = std(eAdj(~hop));
sdCS = sqrt(mean(errCS.^2, 2));  sdLS = sqrt(mean(errLS.^2, 2));
fprintf('rho = %.3f, SNR = %d dB, screen rms = %d ns, Lcorr = %.0f el, %d realizations\n', rho, snrdB, rmsNs, Lcorr, nReal);
fprintf('adjacent-pair jitter sigma_D = %.2f ns (peak hops: %.2f %% of pairs)\n', sigD*1e9, 100*mean(hop));
fprintf('rms error at element 64: cumsum = %.1f ns (random walk sqrt(63)*sigma_D = %.1f ns), weighted LS lags 1..%d = %.1f ns\n', ...
    sdCS(end)*1e9, sqrt(63)*sigD*1e9, nLagLS, sdLS(end)*1e9);
fprintf('phase error at element 64, cumsum: %.0f deg at %.0f MHz\n', sdCS(end)*f0*360, f0/1e6);

%% Figure
figure;
subplot(1, 2, 1);
plot(el, (tau - tau(1))*1e9, 'k-', el, tauCS*1e9, 'o-', el, tauLS*1e9, 's-');
xlabel('element'); ylabel('delay re element 0 (ns)'); grid on;
legend('true screen', 'cumsum', sprintf('WLS, lags 1..%d', nLagLS), 'Location', 'southwest');
title('Last realization');
subplot(1, 2, 2);
plot(el, sdCS*1e9, 'o-', el, sdLS*1e9, 's-', el, sqrt(el)*sigD*1e9, 'k--');
xlabel('element n'); ylabel('rms error (ns)'); grid on;
legend('cumsum', 'weighted LS', 'sqrt(n) \sigma_D', 'Location', 'northwest');
title('Error growth');

What to look for. With \(\rho\) = 0.984 (VCZ value for a uniform 64-element transmit) and 20 dB, \(\sigma_\Delta\) = 1.74 ns, near the decorrelation floor of 1.6 ns. The cumulative-sum error reaches 12.2 ns (22°) at element 64, against \(\sqrt{63}\,\sigma_\Delta\) = 13.8 ns; least squares over lags 1–3 gives 11.7 ns. Set rho = 0.95: lag-3 pairs start to hop and LS becomes worse than the cumulative sum. Set rho = 0.9999, snrdB = 10: noise dominates and LS wins clearly.

m11_block_matching_2d.m2D block matching on RF speckle: separable, axial-first and 2D quadratic subsample fits, NCC vs SSD · change FWHMx, dx, useEnv · runtime ≈ 6 s 
%% m11_block_matching_2d.m
% Purpose : 2D speckle tracking by block matching. An RF speckle image
%           (axial: f0 carrier under a Gaussian envelope; lateral: Gaussian
%           beam, no carrier) is translated by a random subsample shift.
%           A kernel from frame 1 is matched over a search region in
%           frame 2 with 2D NCC and with SSD; the integer peak is refined
%           with (a) separable 1D parabolic fits through the integer peak,
%           (b) axial-first: an axial parabola in each of the 3 lateral
%           columns, then a lateral parabola through the 3 interpolated peak
%           heights, (c) a 2D quadratic least-squares fit (with cross term)
%           on the 3x3 neighbourhood.
% Change  : FWHMx (lateral beam width), dx (line spacing), snrdB,
%           kernel size, useEnv (track the envelope instead of RF).
% Expect  : f0 = 5 MHz, 60 % BW, lateral FWHM 0.6 mm, 20 dB, RF tracking
%           (rms axial / lateral error): separable 0.59 / 15.5 um,
%           axial-first 0.21 / 2.7 um, 2D quadratic LSQ 0.19 / 4.8 um,
%           SSD = NCC. Almost the same at 60 dB SNR: the errors come from
%           the finite speckle kernel and from axial-lateral coupling, not
%           from electronic noise. The separable lateral fit uses the integer
%           axial row, up to half a sample off the RF peak; there the carrier
%           makes the lateral profile of one speckle realization tilted and
%           distorted. Axial-first gives a lateral/axial ratio of about 13,
%           close to the curvature prediction sqrt(R''_ax/R''_lat) = 14.9.
%           Envelope tracking (useEnv = 1) makes the axial error 3-4x worse
%           (no carrier: curvature ratio only 2.6).
rng(1);

%% Parameters
c = 1540;  f0 = 5e6;  fs = 40e6;  B6 = 0.6;    % (m/s), (Hz), (Hz), -6 dB fractional BW
dz = c / (2*fs);                               % axial pixel (m) = 19.25 um
dx = 0.1e-3;                                   % lateral line spacing (m)
FWHMx = 0.6e-3;                                % lateral pulse-echo amplitude FWHM (m)
Nz = 256;  Nx = 64;                            % image size (pixels)
kz = 80;   kx = 15;                            % kernel: 80 axial x 15 lateral pixels (1.54 x 1.5 mm)
sz = 6;    sx = 4;                             % search half-range (pixels)
snrdB = 20;  nTrial = 150;  useEnv = 0;        % SNR, trials, 1 = track |IQ| instead of RF

%% Spectral model (one-sided in axial frequency so that RF and IQ come together)
fz = (0:Nz-1)' / (Nz*dz);  fz(fz >= 1/(2*dz)) = fz(fz >= 1/(2*dz)) - 1/dz;   % axial spatial freq (cycles/m)
fx = (0:Nx-1) / (Nx*dx);   fx(fx >= 1/(2*dx)) = fx(fx >= 1/(2*dx)) - 1/dx;   % lateral spatial freq
fz0 = 2*f0/c;                                  % axial carrier (cycles/m), two-way
sfz = B6*fz0/2/sqrt(2*log(2));                 % axial spectral std (cycles/m)
sigx = FWHMx / (2*sqrt(2*log(2)));             % lateral PSF std (m)
Hz = exp(-(fz - fz0).^2 / (2*sfz^2)) .* (fz > 0);
Hx = exp(-2*pi^2*sigx^2*fx.^2);                % FT of the Gaussian lateral PSF
H  = Hz * Hx;                                  % Nz x Nx
fitPar = @(Rm, R0, Rp) 0.5 * (Rm - Rp) ./ (Rm - 2*R0 + Rp);
% 2D quadratic LSQ on a 3x3 patch: P(u,v) = a + b u + c v + d u^2 + e v^2 + g u v
[uu, vv] = ndgrid(-1:1, -1:1);
Q = [ones(9,1) uu(:) vv(:) uu(:).^2 vv(:).^2 uu(:).*vv(:)];  Qp = (Q'*Q) \ Q';

z0 = 80;  x0 = 25;                             % kernel top-left (1-based)
iz = z0:z0+kz-1;  ix = x0:x0+kx-1;
errSep = zeros(nTrial, 2);  errAF = errSep;  errQ = errSep;  errSSD = errSep;
for t = 1:nTrial
    d = [rand rand];                           % true shift (axial, lateral) in pixels
    W0 = (randn(Nz,Nx) + 1j*randn(Nz,Nx)) .* H;
    a1 = ifft2(W0);
    a2 = ifft2(W0 .* exp(-1j*2*pi*(fz*dz*d(1) + fx*dx*d(2))));
    sc = sqrt(mean(real(a1(:)).^2));
    a1 = a1/sc + (randn(Nz,Nx) + 1j*randn(Nz,Nx)) * 10^(-snrdB/20);   % complex noise: RF part has SNR = snrdB
    a2 = a2/sc + (randn(Nz,Nx) + 1j*randn(Nz,Nx)) * 10^(-snrdB/20);
    if useEnv, I1 = abs(a1); I2 = abs(a2); else, I1 = real(a1); I2 = real(a2); end
    ker = I1(iz, ix);  ker = ker - mean(ker(:));  ek = sum(ker(:).^2);
    R = zeros(2*sz+1, 2*sx+1);  SSD = R;
    for p = -sz:sz
        for q = -sx:sx
            blk = I2(iz+p, ix+q);  b0 = blk - mean(blk(:));
            R(p+sz+1, q+sx+1) = sum(ker(:).*b0(:)) / sqrt(ek * sum(b0(:).^2));
            df = I1(iz, ix) - blk;  SSD(p+sz+1, q+sx+1) = sum(df(:).^2);
        end
    end
    R(:, [1 end]) = -Inf;  R([1 end], :) = -Inf;           % keep the peak off the border
    [~, k] = max(R(:));  [ip, iq] = ind2sub(size(R), k);
    P = R(ip-1:ip+1, iq-1:iq+1);
    sep = [fitPar(P(1,2), P(2,2), P(3,2)), fitPar(P(2,1), P(2,2), P(2,3))];
    % axial-first: axial parabola in each column, lateral parabola through the vertex heights,
    % then the axial offset interpolated to that lateral position (follows a tilted ridge)
    da = fitPar(P(1,:), P(2,:), P(3,:));                     % axial offsets of the 3 columns
    pk = P(2,:) - 0.25 * (P(1,:) - P(3,:)) .* da;            % parabola vertex heights
    u  = fitPar(pk(1), pk(2), pk(3));                        % lateral offset
    af = [da(2) + 0.5*(da(3) - da(1))*u + 0.5*(da(1) - 2*da(2) + da(3))*u^2, u];   % axial offset read at that lateral position
    co = Qp * P(:);                                          % [a b c d e g]
    uv = -[2*co(4) co(6); co(6) 2*co(5)] \ [co(2); co(3)];   % stationary point of the quadratic
    if any(~isfinite(uv)) || any(abs(uv) > 1), uv = sep'; end
    base = [ip - sz - 1, iq - sx - 1];
    errSep(t, :) = base + sep - d;
    errAF(t, :)  = base + af - d;
    errQ(t, :)   = base + uv' - d;
    S2 = -SSD;  S2(:, [1 end]) = -Inf;  S2([1 end], :) = -Inf;
    [~, k] = max(S2(:));  [jp, jq] = ind2sub(size(S2), k);  Ps = -SSD(jp-1:jp+1, jq-1:jq+1);
    errSSD(t, :) = [jp - sz - 1, jq - sx - 1] + [fitPar(Ps(1,2), Ps(2,2), Ps(3,2)), fitPar(Ps(2,1), Ps(2,2), Ps(2,3))] - d;
end

%% Results (um)
pix = [dz dx] * 1e6;
rs = @(e) sqrt(mean(e.^2, 1)) .* pix;        % rms error per direction (um)
modes = {'RF', 'envelope'};
fprintf('pixels: axial %.2f um, lateral %.0f um; lateral FWHM %.2f mm; SNR %d dB; %s tracking\n', pix, FWHMx*1e3, snrdB, modes{useEnv+1});
fprintf('rms error (um)            axial    lateral   ratio\n');
e = rs(errSep); fprintf('NCC + separable parabola  %6.2f   %7.2f   %5.1f\n', e, e(2)/e(1));
e = rs(errAF);  fprintf('NCC + axial-first fits   %6.2f   %7.2f   %5.1f\n', e, e(2)/e(1));
e = rs(errQ);   fprintf('NCC + 2D quadratic LSQ    %6.2f   %7.2f   %5.1f\n', e, e(2)/e(1));
e = rs(errSSD); fprintf('SSD + separable parabola  %6.2f   %7.2f   %5.1f\n', e, e(2)/e(1));
% curvature prediction for RF: -R''(0) axial = (2 pi fz0)^2 + 1/(2 sz^2), lateral = 1/(2 sigx^2)
sigz = 1/(2*pi*sfz);
cAx = (2*pi*fz0)^2 + 1/(2*sigz^2);  cEnv = 1/(2*sigz^2);  cLat = 1/(2*sigx^2);
fprintf('curvature prediction: lateral/axial = %.1f (RF), %.1f (envelope, |IQ| ACF approx.)\n', sqrt(cAx/cLat), sqrt(cEnv/cLat));

%% Figure: last trial
figure('Position', [100 100 1000 420]);
subplot(1, 3, 1);
imagesc(((1:Nx)-1)*dx*1e3, ((1:Nz)-1)*dz*1e3, real(a1)); colormap(gray); hold on;
plot(([x0 x0+kx x0+kx x0 x0]-1)*dx*1e3, ([z0 z0 z0+kz z0+kz z0]-1)*dz*1e3, 'y-', 'LineWidth', 1.5); hold off;
xlabel('lateral (mm)'); ylabel('depth (mm)'); title('RF frame 1 + kernel');
subplot(1, 3, 2);
Rv = R;  Rv(~isfinite(Rv)) = NaN;
imagesc((-sx:sx)*dx*1e3, (-sz:sz)*dz*1e3, Rv); colorbar; hold on;
plot(d(2)*dx*1e3, d(1)*dz*1e3, 'w+', 'MarkerSize', 12, 'LineWidth', 2); hold off;
xlabel('lateral lag (mm)'); ylabel('axial lag (mm)'); title('2D NCC surface');
subplot(1, 3, 3);
plot(errSep(:,2)*pix(2), errSep(:,1)*pix(1), '.', errAF(:,2)*pix(2), errAF(:,1)*pix(1), 'o'); grid on;
xlabel('lateral error (um)'); ylabel('axial error (um)'); title('Error scatter');
legend('separable', 'axial-first', 'Location', 'southoutside');

What to look for. Rms axial / lateral error: separable 0.59 / 15.5 µm, axial-first 0.21 / 2.7 µm, 2D quadratic LSQ 0.19 / 4.8 µm; SSD behaves like NCC. The axial-first lateral/axial ratio (≈ 13) approaches the curvature prediction of 14.9. The errors barely change at 60 dB: they come from the finite kernel and axial–lateral coupling. With useEnv = 1 the axial error grows 3–4×.

11.12Self-check

Q1An ARFI tracker reports a peak displacement of 0.4 samples at \(f_s\) = 40 MHz. What is that in ns of two-way delay and in µm of axial displacement (\(c\) = 1540 m/s)?

0.4 × 25 ns = 10 ns of two-way delay. \(\Delta z = cD/2 = 1540 \times 10\times10^{-9}/2\) = 7.7 µm. Forgetting the factor 1/2 (treating the delay as one-way) doubles the displacement: a common error when porting radar or sonar code.

Q2Compute the Nyquist velocity for \(f_0\) = 5 MHz and PRF = 5 kHz. What is the maximum displacement between pulses before the Kasai phase wraps?

\(v_{\mathrm{Nyq}} = c\cdot\mathrm{PRF}/(4f_0)\) \(= 1540 \cdot 5000 / (4 \cdot 5\times10^6)\) = 0.385 m/s. Per pulse this is \(v_{\mathrm{Nyq}}/\mathrm{PRF}\) = 77 µm = \(\lambda/4\) (\(\lambda\) = 308 µm). In two-way delay that is half an RF period, 100 ns: the phase step reaches \(\pi\).

Q3Why is the bias of every symmetric three-point interpolator exactly zero at \(d\) = 0 and \(d\) = 0.5 samples, and largest in between?

\(R(k) = R_{ss}((k-d)T_s)\) with \(R_{ss}\) even. At \(d = 0\) the samples are symmetric about \(k = 0\) (\(R_- = R_+\)), so any symmetric fit returns 0. At \(d = 0.5\) the samples are symmetric about \(k = 0.5\) (\(R(0) = R(1)\), \(R(-1) = R(2)\)), so the fit returns 0.5. Mirroring the whole configuration about ½ gives \(\mathrm{bias}(d) = -\mathrm{bias}(1-d)\), so the bias is odd about ½. A smooth odd function with zeros at 0, ½ and 1 has its extremes in between: near 0.3 and 0.7 for the RF parabola at \(f_s/f_0\) = 8. (In expectation: one finite speckle kernel is not exactly symmetric.)

Q4For broadband RF speckle at \(f_s/f_0\) = 8, rank parabolic, Gaussian, cosine, upsample ×8 and phase-zero estimators by worst-case bias.

From Demo 11.2 at 5 MHz, 60 %, 40 MHz: phase-zero (0, if \(f_0\) is exact) < upsample ×8 + parabolic (≈ \(4\times10^{-4}\) sample) < cosine (≈ \(1.4\times10^{-3}\)) < parabolic (≈ 0.012) < Gaussian (≈ 0.025). The Gaussian fit models an envelope peak, not an RF cosine peak, so it is the worst choice on RF and the right one on envelope or \(\lvert\mathrm{IQ}\rvert\) correlation. Céspedes et al. (1995) reached the same ordering for parabolic vs cosine on RF echoes. With an \(f_0\) error the phase estimator loses its advantage (5 % error gives about 0.025 sample).

Q5Evaluate the Walker–Trahey bound for \(f_0\) = 5 MHz, \(B\) = 0.6, \(T\) = 2 µs, SNR = 30 dB, \(\rho\) = 1, then for \(\rho\) = 0.99. What does doubling \(T\) do in each case?

\(3/(2 f_0^3 \pi^2 T (B^3+12B))\) = 3/(2 · 1.25·10²⁰ · 9.87 · 2·10⁻⁶ · 7.416) = \(8.2\times10^{-17}\) s². With \(\mathrm{SNR}_p\) = 1000 the bracket is \((1.001)^2 - 1\) = 0.002, so \(\sigma \ge \sqrt{1.64\times10^{-19}}\) = 0.41 ns. With \(\rho\) = 0.99 the bracket is 1.002/0.9801 − 1 = 0.0223, so \(\sigma \ge\) 1.35 ns: decorrelation dominates. Doubling \(T\) divides both by \(\sqrt{2}\) (0.29 ns, 0.96 ns), but only if \(\rho\) does not fall with \(T\), which it does under strain (Section 11.5).

Q6Monte Carlo at 5 dB SNR gives \(\sigma\) = 100 ns while the CRLB is 8 ns. Is the estimator broken?

No. This is the threshold region: some trials pick a side peak one or more periods (200 ns at 5 MHz) away, so the error distribution is a mixture and its variance is dominated by those outliers (e.g. a quarter of trials hopping one period: \(0.25 \times 200^2 \approx (100\ \mathrm{ns})^2\); m11_jitter_vs_crlb.m sees 36 % outliers, some multi-period, and 202 ns at 5 dB). The CRLB only describes local errors. Lengthen the kernel, widen the bandwidth, shrink the search range, use envelope correlation for the coarse search, or reject low-NCC kernels.

Q7Adjacent-element delay estimates have \(\sigma_\Delta\) = 3 ns. After integrating across 64 elements, what is the expected error at the far end, and how would you reduce it?

A random walk: \(\sigma_{64} \approx \sqrt{63} \times 3 \approx\) 24 ns, about 43° of phase at 5 MHz, comparable to a moderate aberrator. Remedies: correlate each element against a beamsum (or a running corrected sum) reference so errors do not accumulate; solve a least-squares problem using lags 1, 2, 3 (overdetermined, which damps the drift); weight each \(\Delta\tau\) by its NCC; iterate the correction, since a better focus raises \(\rho\).

Q8Why does FFT-based correlation of two 256-sample records without zero padding give wrong delays for large shifts?

The inverse DFT of \(X^*Y\) is the circular cross-correlation (Module 04): the lag-\(k\) value includes products of samples that wrapped around from the other end. For large \(k\) most products are wrapped, unrelated content, and the peak can move. Zero-padding to ≥ 511 samples (use 512) makes circular equal linear. Normalizing by \(\sqrt{\text{energy}}\) of the actually overlapping segments gives true NCC.

11.13Figures from the literature

Measured or benchmark data to hold this module's models against. Images are not reproduced; open the linked figure. Figure numbers were checked against the linked full text or the publisher's figure list.

Inter-element coherence versus lag: theory against in vivo tissue

M. A. Lediju, G. E. Trahey, B. C. Byram, J. J. Dahl, "Short-lag spatial coherence of backscattered echoes: imaging characteristics," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58(7):1377–1388, 2011. doi:10.1109/TUFFC.2011.1957 (author copy: PULSE Lab). Fig. 1.

Fig. 1 plots the correlation of RF echoes between receive elements against their separation (lag) for a point target, for diffuse speckle and for in vivo thyroid tissue. Diffuse scatter gives the van Cittert–Zernike triangle falling from 1 at zero lag to 0 at lag \(N-1\); the thyroid curve falls below it, which the authors attribute to clutter, phase aberration and noise, with the largest coherence losses at short lags. That short-lag value is exactly the adjacent-element \(\rho\) that sets the decorrelation floor of (11.13) in phase-screen estimation: the theoretical \(1 - 1/N\) is a best case that real tissue does not reach.

Subpixel localization errors of seven ULM algorithms, benchmarked

B. Heiles, A. Chavignon, V. Hingot et al., "Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy," Nature Biomedical Engineering 6:605–616, 2022. doi:10.1038/s41551-021-00824-8. Fig. 4 and Fig. 5.

Fig. 4 shows the directional localization errors and statistical indices of each algorithm on the "in silico flow" dataset at 30 dB; Fig. 5 studies how a non-uniform PSF in different directions affects localization. Read them as the 2D, single-frame version of Section 11.2: a localizer is a subsample peak estimator on an envelope image, its errors depend on the assumed peak shape (weighted average, interpolation, Gaussian fit), and directional bias toward pixel centres is the localization analogue of pixel locking.

Localization and frame-to-frame tracking of bubbles in the rat brain

C. Errico, J. Pierre, S. Pezet et al., "Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging," Nature 527:499–502, 2015. doi:10.1038/nature16066. Fig. 1 and Fig. 2.

Fig. 1 gives the principle of ultrafast ULM (plane-wave acquisition, clutter rejection, localization, tracking); Fig. 2 quantifies spatial resolution and flow in the rat cortex through a thinned skull. The pairing step of Section 11.6(d) is what turns localizations into the velocities of Fig. 2: at kilohertz frame rates bubbles move a small fraction of a wavelength per frame, so the search radius, not the localization precision, decides whether neighbours are mispaired.

Phase-shift tracking between transmit angles for sound-speed maps

P. Stähli, M. Kuriakose, M. Frenz, M. Jaeger, "Forward model for quantitative pulse-echo speed-of-sound imaging," arXiv:1902.10639v2 (2019); published as "Improved forward model for quantitative pulse-echo speed-of-sound imaging," Ultrasonics 108, 2020. Fig. 2 and Fig. 7 (arXiv numbering).

The delay estimator here is the phase of the kernel-averaged Hermitian product of complex images from two transmit angles (2 × 2 mm kernel), tracked in 2° steps and accumulated over 10° to stay inside the \(\pm\pi\) wrap; Fig. 2 relates the measured phase shift to the shift of the echo-envelope centroid, and Fig. 7 shows the resulting B-mode and speed-of-sound images of a volunteer's abdominal wall and liver. Every design choice maps onto this module: kernel size against the bound (11.13), angle step against the wrap (11.16), and inter-angle decorrelation as the \(\rho\) that limits the jitter.

11.14Go deeper

  • S. M. Kay, Fundamentals of Statistical Signal Processing, Vol. I: Estimation Theory, Prentice Hall, 1993: Ch. 3 (CRLB), Ch. 7 (MLE). Batch 2's estimation module builds on these.
  • W. F. Walker, G. E. Trahey, "A fundamental limit on delay estimation using partially correlated speckle signals," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 42(2), 301–308, 1995.
  • C. H. Knapp, G. C. Carter, "The generalized correlation method for estimation of time delay," IEEE Trans. Acoust. Speech Signal Process., 24(4), 320–327, 1976.
  • G. C. Carter, "Coherence and time delay estimation," Proc. IEEE, 75(2), 236–255, 1987.
  • A. J. Weiss, E. Weinstein, "Fundamental limitations in passive time delay estimation, Part I: Narrow-band systems," IEEE Trans. Acoust. Speech Signal Process., 31(2), 472–486, 1983 (Ziv–Zakai / Barankin and the threshold).
  • J. Ziv, M. Zakai, "Some lower bounds on signal parameter estimation," IEEE Trans. Inf. Theory, 15(3), 386–391, 1969.
  • I. Céspedes, Y. Huang, J. Ophir, S. Spratt, "Methods for estimation of subsample time delays of digitized echo signals," Ultrasonic Imaging, 17(2), 142–171, 1995.
  • F. Viola, W. F. Walker, "A comparison of the performance of time-delay estimators in medical ultrasound," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 50(4), 392–401, 2003.
  • G. F. Pinton, J. J. Dahl, G. E. Trahey, "Rapid tracking of small displacements with ultrasound," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 53(6), 1103–1117, 2006.
  • C. Kasai, K. Namekawa, A. Koyano, R. Omoto, "Real-time two-dimensional blood flow imaging using an autocorrelation technique," IEEE Trans. Sonics Ultrason., 32(3), 458–464, 1985.
  • T. Loupas, J. T. Powers, R. W. Gill, "An axial velocity estimator for ultrasound blood flow imaging, based on a full evaluation of the Doppler equation by means of a two-dimensional autocorrelation approach," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 42(4), 672–688, 1995.
  • T. Varghese, J. Ophir, "A theoretical framework for performance characterization of elastography: the strain filter," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 44(1), 164–172, 1997.
  • S. W. Flax, M. O'Donnell, "Phase-aberration correction using signals from point reflectors and diffuse scatterers: basic principles," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 35(6), 758–767, 1988.
  • M. Jaeger, G. Held, S. Peeters, S. Preisser, M. Grünig, M. Frenz, "Computed ultrasound tomography in echo mode for imaging speed of sound using pulse-echo sonography: proof of principle," Ultrasound Med. Biol., 41(1), 235–250, 2015.
  • C. Errico et al., "Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging," Nature, 527, 499–502, 2015.