Part II · Module 12

Estimation theory and the Cramér–Rao bound

Every number your pipeline produces is an estimate: a bubble position in ULM, an element-to-element delay, a Doppler frequency, a sound speed, an attenuation slope. This module builds the language for how precise those numbers can be (likelihood, Fisher information, the Cramér–Rao bound) and for the three ways the bound misleads: thresholds, nuisance parameters and biased estimators.

You will be able to

Split any error into bias and variance, and say which estimators are efficient for a given noise model
Write the Fisher information for a signal in real or complex Gaussian noise, and get the factor 2 right
Derive \(\sigma_x \ge w/\sqrt{\mathrm{ENR}}\) for ULM localization, convert it to µm and λ, and judge a localizer against it
Recognize the threshold region, and quantify what unknown nuisance parameters cost (neighbour bubbles, \(t_0\) in sound-speed estimation)
Explain when a biased or Bayesian estimator legitimately beats the CRLB in MSE

Assumes

Module 05 · Matrices as operators: least squares as a projection
Module 07 · Random vectors, covariance & correlation: circular complex Gaussian noise and speckle
Module 08 · Analytic signal & IQ demodulation: IQ data, phase of a complex signal
Module 11 · Correlation & time-delay estimation: delay estimation, Walker–Trahey bound, false peaks

Exact formulaclosed form plotted as isSimulationMonte Carlo or matrix solve run here (JS in the page, or the named Octave script)Approx. modelsimplified physics (Gaussian PSF, white circular noise, straight-ray delays); assumptions stated under the badgeSchematicdrawing for intuition, not computedNo curve on this page is hand-drawn: every plotted line comes from the formula or the simulation named next to it.

12.1Estimators, bias, variance and MSE

An estimator \(\hat\theta = g(\mathbf{y})\) is a function of noisy data, so it is a random variable. Three numbers describe it: its bias \(b(\theta) = \mathrm{E}[\hat\theta] - \theta\), its variance \(\operatorname{var}(\hat\theta) = \mathrm{E}[(\hat\theta - \mathrm{E}\hat\theta)^2]\), and their combination, the mean squared error. Adding and subtracting \(\mathrm{E}\hat\theta\) inside the square, the cross term has zero mean:

\[ \mathrm{MSE}(\hat\theta) = \mathrm{E}\big[(\hat\theta - \theta)^2\big] = b^2(\theta) + \operatorname{var}(\hat\theta) \]

(12.1)

\(\mathrm{MSE}\): mean-squared error of the estimator, units of \(\theta^2\)
\(\theta,\ \hat\theta\): true parameter and its estimate
\(\mathrm{E}[\cdot]\): expectation over the noise
\(b(\theta)\): bias \(\mathrm{E}[\hat\theta] - \theta\)
\(\operatorname{var}(\hat\theta)\): estimator variance

Bias and variance trade against each other. An RMSE of 3.6 µm can be 2 µm of bias plus 3 µm of jitter: \(\sqrt{2^2 + 3^2} = 3.6\) µm. Report both parts, because only the variance averages down with more data.

An estimator is consistent if \(\hat\theta \to \theta\) as \(N \to \infty\), unbiased if \(b \equiv 0\) for every \(\theta\), and efficient if it is unbiased and its variance equals the Cramér–Rao lower bound (Section 12.3). The bound is a property of the data model, not of any estimator: change the noise law and the best estimator changes.

The simplest case shows this. Estimate a constant \(A\) from \(N\) samples \(y_n = A + w_n\), noise variance \(\sigma^2\):

Noise law (variance \(\sigma^2\))	CRLB on \(\operatorname{var}(\hat A)\)	Sample mean	Sample median	Midrange (max+min)/2
Gaussian	\(\sigma^2/N\)	\(\sigma^2/N\) (efficient)	\(\approx \pi\sigma^2/(2N)\) (64 % efficiency)	poor, \(\propto 1/\ln N\)
Laplace	\(\sigma^2/(2N)\)	\(\sigma^2/N\) (50 %)	\(\to \sigma^2/(2N)\) (asymptotically efficient)	does not decrease with N
Uniform	does not exist	\(\sigma^2/N\)	\(3\sigma^2/N\)	\(6\sigma^2/((N+1)(N+2)) \propto 1/N^2\)

The uniform row breaks the regularity condition: the support of \(p(y; A)\) moves with \(A\), so the derivative of the log-likelihood is not defined at the edges and the bound does not apply. The midrange then beats \(1/N\). Nothing in ultrasound IQ noise behaves like that, but the lesson carries: a bound is only as good as its assumptions.

Demo 12.1

Mean, median and midrange against the CRLB

SimulationMonte Carlo in the page (JS, fixed seed): 4000 trials per N (2000 for N > 60); CRLB curve from the closed form.

sample meansample medianmidrangeCRLB density N(0, N·CRLB/σ²)

Noise law (unit variance)

Samples N25

Try: top panel shows \(\sqrt{N}(\hat A - A)/\sigma\) from 4000 trials (2000 for N > 60), so an efficient estimator always matches the black curve. Gaussian: the mean sits on the bound and the median is wider by \(\sqrt{\pi/2} = 1.25\) in std. Switch to Laplace: the median approaches the bound \(\sigma^2/(2N)\) slowly (\(N\cdot\operatorname{var} \approx 0.68\) at N = 25, 0.6 at N = 129) while the mean stays at 1. Uniform: no bound exists, and the midrange curve collapses as N grows because its variance falls as \(1/N^2\).

Ultrasound link

Speckle intensity \(I\) is exponential with mean \(\mu\) (Module 07). The sample mean of intensities is the MLE of \(\mu\) and is efficient (\(\operatorname{var} = \mu^2/N\)). Averaging log-compressed values instead estimates \(\ln\mu\) with a bias of \(-\gamma\) (−2.51 dB in \(10\log_{10}\)) and a variance \(\pi^2/6 = 1.64\) times the CRLB: 61 % efficiency. Backscatter and attenuation estimators should average intensity (or power spectra), not decibels.

12.2Likelihood, the MLE and least squares

The likelihood is the data density read as a function of the parameter: \(L(\theta) = p(\mathbf{y}; \theta)\). The maximum-likelihood estimator is \(\hat\theta_{\mathrm{ML}} = \arg\max_\theta L(\theta)\). Under regularity conditions it has four properties that make it the default:

Invariance: the MLE of \(\alpha = g(\theta)\) is \(g(\hat\theta_{\mathrm{ML}})\), for any \(g\).
Consistency, and asymptotic normality: \(\hat\theta_{\mathrm{ML}} \sim \mathcal{N}(\theta, I^{-1}(\theta))\) as \(N \to \infty\) or as SNR → ∞.
Asymptotic efficiency: it reaches the CRLB in that limit.
If an efficient estimator exists at finite N, the MLE is it.

None of these promise unbiasedness at finite N or low SNR. Invariance is the classic trap: if \(\hat a\) is the MLE of a complex amplitude, \(\lvert\hat a\rvert^2\) is the MLE of the power, and it is biased upward by the noise term, \(\mathrm{E}\lvert\hat a\rvert^2 = \lvert a\rvert^2 + \sigma^2/\lVert\mathbf{s}\rVert^2\).

Least squares is the MLE in Gaussian noise

For real white Gaussian noise of variance \(\sigma^2\), and for circular complex Gaussian noise with \(\mathrm{E}\lvert w\rvert^2 = \sigma^2\):

\[ \begin{aligned} &\text{real:} \\ &\ln p = -\frac{1}{2\sigma^2}\,\lVert \mathbf{y} - \mathbf{s}(\theta) \rVert^2 + \text{const} \\[4pt] &\text{complex:} \\ &\ln p = -\frac{1}{\sigma^2}\,\lVert \mathbf{y} - \mathbf{s}(\theta) \rVert^2 - N \ln(\pi\sigma^2) \end{aligned} \]

(12.2)

\(\ln p\): log-likelihood of the data given \(\theta\)
\(\mathbf{y}\): measured data vector (\(N\) samples)
\(\mathbf{s}(\theta)\): noise-free signal model for parameter \(\theta\)
\(\sigma^2\): noise variance per real sample (real case) or total per complex sample (complex case)
\(N\): number of samples
\(\lVert\cdot\rVert\): Euclidean norm
\(\mathrm{const}\): terms independent of \(\theta\)

Maximizing either is minimizing the squared residual: (nonlinear) least squares. With coloured noise of covariance \(\mathbf{C}\), the exponent becomes \((\mathbf{y} - \mathbf{s})^{H}\mathbf{C}^{-1}(\mathbf{y} - \mathbf{s})\): weighted least squares, or pre-whiten with \(\mathbf{C}^{-1/2}\) and then use ordinary least squares.

For the linear model \(\mathbf{y} = \mathbf{H}\theta + \mathbf{w}\) the MLE is the least-squares solution \(\hat\theta = (\mathbf{H}^{T}\mathbf{H})^{-1}\mathbf{H}^{T}\mathbf{y}\) (Module 05), with covariance \(\sigma^2(\mathbf{H}^{T}\mathbf{H})^{-1}\). That covariance equals the CRLB exactly, at any N: in the linear Gaussian model least squares is efficient, not just asymptotically.

Unknown complex amplitude: the matched filter scans the parameter

Most ultrasound models have the form \(\mathbf{s} = a\,\mathbf{g}(\theta)\) with a complex amplitude \(a\) (echo strength and phase) that nobody cares about. Minimizing \(\lVert\mathbf{y} - a\mathbf{g}\rVert^2\) over \(a\) gives \(\hat a = \mathbf{g}^{H}\mathbf{y}/\lVert\mathbf{g}\rVert^2\), and substituting back leaves a function of \(\theta\) alone:

\[ \hat\theta_{\mathrm{ML}} = \arg\max_{\theta}\ \frac{\lvert \mathbf{g}(\theta)^{H}\mathbf{y} \rvert^2}{\lVert \mathbf{g}(\theta) \rVert^2} \]

(12.3)

\(\hat\theta_{\mathrm{ML}}\): maximum-likelihood estimate of \(\theta\)
\(\mathbf{g}(\theta)\): signal shape for parameter \(\theta\) (tone, PSF centred at \(x_0\), delayed pulse)
\(\mathbf{y}\): measured data vector
\((\cdot)^{H}\): conjugate transpose
\(\lvert\cdot\rvert,\ \lVert\cdot\rVert\): magnitude and Euclidean norm
\(\arg\max_\theta\): value of \(\theta\) that maximizes the expression
\(\hat a\): profiled-out complex amplitude, \(\mathbf{g}^{H}\mathbf{y}/\lVert\mathbf{g}\rVert^2\)

The concentrated (profile) likelihood. For a tone, \(\mathbf{g}(\omega) = e^{j\omega n}\) and this is the periodogram: its peak is the frequency MLE. For a bubble, \(\mathbf{g}(x_0)\) is the PSF centred at \(x_0\): the MLE slides a matched filter across the image. For a delay, it is the cross-correlation peak (Module 11).

12.3Fisher information and the Cramér–Rao bound

The score \(\partial \ln p(\mathbf{y}; \theta)/\partial\theta\) has zero mean when the support of \(p\) does not depend on \(\theta\). Its variance is the Fisher information, which also equals the expected curvature of the log-likelihood:

\[ \begin{aligned} I(\theta) &= \mathrm{E}\!\left[\left(\frac{\partial \ln p}{\partial\theta}\right)^{2}\right] = -\mathrm{E}\!\left[\frac{\partial^2 \ln p}{\partial\theta^2}\right] \\ \operatorname{var}(\hat\theta) &\ge \frac{1}{I(\theta)} \quad \text{(unbiased)} \\ \operatorname{var}(\hat\theta) &\ge \frac{\big(1 + b'(\theta)\big)^2}{I(\theta)} \quad \text{(bias } b\text{)} \end{aligned} \]

(12.4)

\(p = p(\mathbf{y}; \theta)\): data density (likelihood)
\(I(\theta)\): Fisher information, 1/units² of \(\theta\)
\(\hat\theta\): estimator of \(\theta\)
\(\operatorname{var}(\hat\theta)\): estimator variance, units² of \(\theta\)
\(\mathrm{E}[\cdot]\): expectation over the data
\(b(\theta),\ b'(\theta)\): estimator bias and its slope with respect to the true value (dimensionless)
\(\partial\): partial derivative

Equality holds if and only if the score factors as \(I(\theta)\,(g(\mathbf{y}) - \theta)\); then \(g\) is the efficient estimator. Information from independent samples adds.

With several parameters the information becomes a matrix and the bound a matrix inequality. The bound on one parameter is a diagonal entry of the inverse, and a function \(\alpha = f(\theta)\) inherits the bound through its gradient:

\[ \begin{aligned} [\mathbf{I}(\theta)]_{ij} &= -\mathrm{E}\!\left[\frac{\partial^2 \ln p}{\partial\theta_i\,\partial\theta_j}\right] \\ \operatorname{cov}(\hat\theta) - \mathbf{I}^{-1} &\succeq 0 \\ \Rightarrow\ \operatorname{var}(\hat\theta_i) &\ge [\mathbf{I}^{-1}]_{ii} \\ \operatorname{var}(\hat\alpha) &\ge \nabla f^{T}\,\mathbf{I}^{-1}\,\nabla f \end{aligned} \]

(12.5)

\(\mathbf{I}(\theta),\ I_{ij}\): Fisher information matrix and its \((i,j)\) entry, as in (12.4)
\(\theta_i,\ \theta_j\): entries of the parameter vector \(\theta\)
\(\operatorname{cov}(\hat\theta)\): covariance matrix of the estimator
\(\succeq 0\): positive semidefinite
\(\operatorname{var}(\hat\theta_i)\): variance of the \(i\)th estimate
\([\mathbf{I}^{-1}]_{ii}\): \(i\)th diagonal entry of the inverse Fisher matrix
\(\alpha = f(\theta)\): function of the parameters, with estimate \(\hat\alpha\)
\(\nabla f\): gradient of \(f\) with respect to \(\theta\)
\((\cdot)^{T}\): transpose

The Gaussian signal-in-noise formula, real and complex

For a deterministic signal \(\mathbf{s}(\theta)\) in Gaussian noise the expectation can be done once and for all:

\[ \begin{aligned} &\text{real noise, variance } \sigma^2: \\ &I_{ij} = \frac{1}{\sigma^2}\sum_n \frac{\partial s[n]}{\partial\theta_i}\,\frac{\partial s[n]}{\partial\theta_j} \\[4pt] &\text{coloured, covariance } \mathbf{C}: \\ &\mathbf{I} = \mathbf{J}^{T}\mathbf{C}^{-1}\mathbf{J} \end{aligned} \]

(12.6)

\(I_{ij},\ \mathbf{I}\): Fisher information entry and matrix, as in (12.5)
\(\sigma^2\): variance of the real white noise per sample
\(s[n]\): noise-free signal sample \(n\), depends on \(\theta\)
\(\theta_i,\ \theta_j\): parameters, as in (12.5)
\(n\): sample index
\(\mathbf{J}\): Jacobian \(\partial\mathbf{s}/\partial\theta\), \(N \times P\) (samples × parameters)
\(\mathbf{C}\): noise covariance, \(N \times N\)
\((\cdot)^{T}\): transpose

\[ \begin{aligned} &\text{circular complex noise, } \mathrm{E}\lvert w\rvert^2 = \sigma^2: \\ &I_{ij} = \frac{2}{\sigma^2}\,\mathrm{Re}\sum_n \left(\frac{\partial s[n]}{\partial\theta_i}\right)^{\!*}\frac{\partial s[n]}{\partial\theta_j} \\[4pt] &\text{coloured, covariance } \mathbf{C}: \\ &\mathbf{I} = 2\,\mathrm{Re}\{\mathbf{J}^{H}\mathbf{C}^{-1}\mathbf{J}\} \end{aligned} \]

(12.7)

\(I_{ij},\ \mathbf{I}\): Fisher information entry and matrix, as in (12.5)
\(\sigma^2\): total noise power per complex sample (\(\sigma^2/2\) per quadrature)
\(w\): circular complex noise sample
\(s[n]\): noise-free complex signal sample
\(\theta_i,\ \theta_j\): parameters
\(n\): sample index
\(\mathbf{J},\ \mathbf{C}\): Jacobian and noise covariance, as in (12.6)
\(\mathrm{Re}\): real part
\((\cdot)^{*},\ (\cdot)^{H}\): conjugate and conjugate transpose
\(\mathrm{E}[\cdot]\): expectation

Information is the energy of the signal's derivative, measured against the noise: a parameter is well determined when changing it changes the signal a lot compared with the noise. The bound depends on the noise only through its covariance, and the parameters only through \(\mathbf{J}\).

Where the factor 2 comes from. A complex sample is two real samples. With total noise power \(\sigma^2 = \mathrm{E}\lvert w\rvert^2\), each quadrature has variance \(\sigma^2/2\), so \(I = \sum(\partial\,\mathrm{Re}\,s)^2/(\sigma^2/2)\) \(+ \sum(\partial\,\mathrm{Im}\,s)^2/(\sigma^2/2)\) \(= (2/\sigma^2)\sum\lvert\partial s\rvert^2\). Three consequences, all checked in m12_crlb_amplitude_mle.m:

Real amplitude of a known complex pulse (phase known), complex noise: \(\operatorname{var}(\hat a) \ge \sigma^2/(2\lVert\mathbf{s}\rVert^2)\).
Complex amplitude (phase unknown): \(\mathrm{E}\lvert\hat a - a\rvert^2 \ge \sigma^2/\lVert\mathbf{s}\rVert^2\), i.e. \(\sigma^2/(2\lVert\mathbf{s}\rVert^2)\) per real component.
Real RF pulse in real noise of variance \(\sigma^2\): \(\operatorname{var}(\hat a) \ge \sigma^2/\lVert\mathbf{s}\rVert^2\).

Demodulation to IQ is an invertible map of the band-limited data, so it neither creates nor destroys information. Factor-2 errors arise when one side of a comparison counts \(\sigma^2\) per quadrature and the other counts it per complex sample. Write down which convention every SNR uses.

Pitfall

The CRLB bounds unbiased estimators only, and only locally. An estimator that clips its output to a search range, or a localizer that returns the pixel centre when the fit fails, is biased; its variance can sit below \(1/I\) and that proves nothing. Always plot bias and RMSE separately against the true parameter.

12.4Frequency estimation, the threshold and global bounds

Slow-time Doppler is a frequency-estimation problem. Model one complex tone in circular complex Gaussian noise, \(y[n] = A\,e^{j(\omega n + \varphi)} + w[n]\), \(n = 0,\dots,N-1\), \(\mathrm{SNR} = A^2/\sigma^2\). With \(\omega\) and \(\varphi\) unknown, the Fisher matrix and its inverse (Schur complement on \(\varphi\)) are:

\[ \begin{aligned} \mathbf{I} &= \frac{2A^2}{\sigma^2}\begin{bmatrix} \sum n^2 & \sum n \\ \sum n & N \end{bmatrix} \\ \Rightarrow\ \operatorname{var}(\hat\omega) &\ge \frac{\sigma^2}{2A^2}\,\frac{N}{N\sum n^2 - \left(\sum n\right)^2} \\ &= \frac{6}{\mathrm{SNR}\cdot N(N^2 - 1)} \\ &\phantom{=}\ \ (\text{rad/sample})^2 \end{aligned} \]

(12.8)

\(\mathbf{I}\): Fisher matrix for the unknowns \((\omega, \varphi)\)
\(A,\ \varphi\): tone amplitude and unknown phase, rad
\(\omega,\ \hat\omega\): true and estimated frequency, rad/sample (Doppler: \(2\pi f_D/\mathrm{PRF}\))
\(n\): pulse index, \(0,\dots,N-1\)
\(N\): ensemble length (pulses)
\(\sigma^2\): total complex noise power per sample
\(\mathrm{SNR}\): per-sample \(A^2/\sigma^2\), linear
\(\sum\): sum over \(n = 0,\dots,N-1\)

\(N\sum n^2 - (\sum n)^2 = N^2(N^2 - 1)/12\) for \(n = 0,\dots,N-1\). The variance falls as \(N^{-3}\): N samples of SNR each, times an \(N^2\) lever arm, because a frequency error produces a phase error that grows linearly along the record.

For a real sinusoid \(x[n] = A\cos(2\pi f n + \varphi) + w[n]\) in real noise of variance \(\sigma^2\) (Kay, Vol. I, Example 3.14), with \(\eta = A^2/(2\sigma^2)\):

\[ \begin{aligned} \operatorname{var}(\hat f) &\ge \frac{12}{(2\pi)^2\,\eta\,N(N^2 - 1)} \\ &\phantom{\ge}\ \ (\text{cycles/sample})^2 \\ \text{i.e.}\quad \operatorname{var}(\hat\omega) &\ge \frac{12}{\eta\,N(N^2 - 1)} \end{aligned} \]

(12.9)

\(\hat f,\ \hat\omega\): estimated frequency (cycles/sample) and angular frequency (rad/sample), \(\omega = 2\pi f\)
\(\eta\): real-sinusoid SNR \(A^2/(2\sigma^2)\), linear
\(N\): number of samples
\(\operatorname{var}\): estimator variance

In each case "SNR" is signal power over noise power. At equal SNR the complex tone has half the variance: the factor 2 again. Kay's form is approximate (it neglects the interaction with the negative-frequency image) and holds for \(f\) not near 0 or ½.

Converting to Doppler: \(f_D = \omega\cdot\mathrm{PRF}/(2\pi)\) and \(v = c f_D/(2f_0)\), so \(\sigma_v = c\cdot\mathrm{PRF}\cdot\sigma_\omega/(4\pi f_0)\). With M = 16 pulses at 10 dB per pulse, \(\sigma_\omega \ge 0.0121\) rad/sample; at PRF = 5 kHz and \(f_0\) = 5 MHz that is 9.65 Hz, or 1.49 mm/s.

For the complex tone the MLE is exactly the periodogram peak (the concentrated likelihood of Section 12.2). The lag-one autocorrelation (Kasai, Module 11) is not efficient: at high SNR its phase estimate averages \(N - 1\) phase differences, which telescope to \((\varphi_{N-1} - \varphi_0)/(N - 1)\), so \(\operatorname{var} \to 1/(\mathrm{SNR}\,(N - 1)^2)\), a factor \(N(N + 1)/(6(N - 1))\) above the bound (5.7 at N = 32). It wins on cost and on robustness when the Doppler spectrum is broad.

The threshold effect and global bounds

The CRLB is local: it describes small errors around the true peak. When noise lifts a periodogram bin elsewhere above the signal bin, the estimate jumps to a random frequency, and the MSE rises toward the variance of a uniform guess, \(\pi^2/3\). The MSE curve has three regions: the asymptotic region on the bound, a threshold region where outliers start, and a no-information region. For the periodogram the outlier probability depends on the integrated SNR \(N\cdot\mathrm{SNR}\) (Rife and Boorstyn, 1974): in m12_frequency_crlb.m (1000 trials), N = 32 sits on the bound for \(N\cdot\mathrm{SNR} \ge 13\) dB and is about 180 times above it at 11 dB, where 2 % of trials are outliers. The 160-trial demo below finds the same knee, but at 11 dB it sees only a handful of outliers, so its ratio there swings by an order of magnitude from seed to seed (14× with the default seed). Global bounds that include these ambiguity errors (Ziv–Zakai, Barankin, Weiss–Weinstein) reproduce the knee; the CRLB never does. Module 11 showed the same mechanism for delays, where the false peaks sit one RF period apart.

Demo 12.2

Periodogram MLE and lag-one estimator against the frequency CRLB

SimulationMonte Carlo in the page (JS), 160 trials per SNR; bound from (12.8).

periodogram peak (MLE)lag-one (Kasai)CRLB 6/(SNR·N(N²−1))Kasai high-SNR limituniform guess π/√3

Ensemble length N

PRF (for velocity)5.0 kHz

Try: \(f_0\) = 5 MHz, c = 1540 m/s, true \(\omega\) uniform over \((-\pi, \pi]\). MLE = zero-padded periodogram peak refined on a fine grid. At N = 32 the MLE leaves the bound between 0 and −4 dB per sample (13 → 11 dB integrated): the knee is the threshold. Raise N to 128: the threshold moves to lower per-sample SNR but stays near the same integrated \(N\cdot\mathrm{SNR}\) (it creeps up slowly with N because there are more bins that can win). The lag-one estimator approaches its high-SNR limit \(N(N+1)/(6(N-1))\) times the CRLB only above about 20 dB.

12.5ULM localization: the bound, the estimators, the pixel

An isolated microbubble in a beamformed IQ image is a complex amplitude times the PSF, plus noise. Near the focus, model the PSF magnitude as separable Gaussians with amplitude standard deviations \(w_x\), \(w_z\) (\(w = \mathrm{FWHM}/2.355\), FWHM of |IQ|):

\[ \begin{aligned} y(x, z) &= a\, g_x(x - x_0)\, g_z(z - z_0) + n(x, z) \\ g_x(u) &= \exp\!\left(-\frac{u^2}{2w_x^2}\right), \quad a \in \mathbb{C}\ \text{unknown} \end{aligned} \]

(12.10)

\(y(x,z)\): complex image around the bubble
\(x,\ z\): lateral and axial position, µm
\(x_0,\ z_0\): true bubble position (the unknowns), µm
\(a\): unknown complex amplitude
\(g_x,\ g_z\): lateral and axial Gaussian PSF profiles
\(u\): dummy position argument, µm
\(w_x,\ w_z\): lateral and axial Gaussian std of the PSF (\(\mathrm{FWHM}/2.355\)), µm
\(n(x,z)\): circular complex Gaussian noise, density \(N_0\) per unit area
\(\mathbb{C}\): complex numbers

Beamformed noise is band-limited, not independent from pixel to pixel; the bound below needs only that the noise is white where the signal has energy.

The derivative with respect to \(x_0\) is \(-a\,g_x' g_z\). For a Gaussian, \(\int g_x'^2\,dx \big/ \int g_x^2\,dx = 1/(2w_x^2)\), and every cross term (with \(\mathrm{Re}\,a\), \(\mathrm{Im}\,a\), \(z_0\)) integrates to zero by symmetry, so the Fisher matrix is diagonal:

\[ \begin{aligned} I_{x_0x_0} &= \frac{2\lvert a\rvert^2}{N_0}\iint \left\lvert\frac{\partial g}{\partial x}\right\rvert^2 \\ &= \frac{2\lvert a\rvert^2}{N_0}\cdot\frac{\pi w_x w_z}{2w_x^2} = \frac{\mathrm{ENR}}{w_x^2} \\ \Rightarrow\ \sigma_x &\ge \frac{w_x}{\sqrt{\mathrm{ENR}}} = \frac{\mathrm{FWHM}_x}{2.355\sqrt{\mathrm{ENR}}} \\ \sigma_z &\ge \frac{w_z}{\sqrt{\mathrm{ENR}}} \end{aligned} \]

(12.11)

\(I_{x_0x_0}\): Fisher information for the lateral position, 1/µm²
\(a,\ g,\ w_x,\ w_z\): amplitude and PSF \(g = g_x g_z\), as in (12.10)
\(N_0\): noise power density per unit image area
\(\iint\): integral over the image plane
\(\partial g/\partial x\): lateral derivative of the PSF
\(\mathrm{ENR}\): energy-to-noise ratio \(\lvert a\rvert^2\pi w_x w_z/N_0\): matched-filter output SNR (linear, not dB)
\(\sigma_x,\ \sigma_z\): lower bounds on the position standard deviation, µm
\(\mathrm{FWHM}_x\): lateral full width at half maximum of the PSF, µm

For noise white over a rectangular band whose edges sit where each 1D PSF energy spectrum falls to −40 dB (\(\lvert k\rvert \le 0.483/w\) cycles per unit length on each axis), the peak-pixel SNR is \(10\log_{10}(\ln 10^4/\pi)\) = 4.7 dB below the ENR in the continuous limit (4.8–5 dB on the discrete grids used here). The offset depends on the band shape: the elliptical −40 dB contour of the 2D spectrum gives 3.6 dB. In real noise with \(\mathrm{ENR}_r = \sum s^2/\sigma^2\), the same steps give \(\sigma_x \ge \sqrt{2}\,w_x/\sqrt{\mathrm{ENR}_r}\).

ENR	\(\sigma_x/\mathrm{FWHM}_x\)	FWHM = λ at 5 MHz (λ = 308 µm)	FWHM = λ at 15 MHz (λ = 103 µm)
12.6 dB	1/10	30.8 µm	10.3 µm
20 dB	1/23.5	13.1 µm	4.4 µm
30 dB	1/74.5	4.1 µm	1.4 µm
40 dB	1/235	1.3 µm	0.44 µm

SchematicDrawing for intuition. Grid, PSF footprint and dimensions to scale for 5 MHz (λ = 308 µm); the σ ellipse is enlarged ×4; pixel shading and profiles use the Gaussian PSF of the text.

Localization geometry. The beamformed bubble echo covers a PSF footprint (solid: −6 dB, dashed: −20 dB) sampled by λ/4 pixels. The MLE slides the template g(x − θ) over the data; the information sits on the flanks of the PSF, where ∂g/∂x is large, so the precision is set by w and the ENR, not by the pixel.

Precision is not resolution. A λ/10 figure in a ULM paper is the standard deviation of one isolated bubble's position. Two bubbles inside one PSF at the same time are a different estimation problem with more unknowns (Section 12.6), and their positions degrade sharply below about one FWHM of separation. Final ULM maps add tracking errors, motion, SVD-filter residue and the finite number of localizations per pixel on top of this bound.

The pixel does not set the bound. \(\sigma_x\) depends on the PSF width and on the noise density. Pixels finer than the Nyquist spacing of the noise band add samples whose noise is correlated, not new information. What the pixel size changes is what a given estimator can extract. A simulation that adds independent noise to every pixel at a fixed per-pixel SNR makes ENR grow as \(1/p^2\), so the "bound" improves with finer pixels: an artefact of the noise model.

How the common localizers compare

Maximum likelihood: maximize \(\lvert\langle y,\, g(\cdot - x_0, \cdot - z_0)\rangle\rvert^2\) over \((x_0, z_0)\), i.e. a Gaussian fit with complex amplitude (12.3). Efficient above threshold at any pixel size that samples the PSF.
Intensity-weighted centroid of \(\lvert y\rvert^2\) in a box of ±1 FWHM: at high ENR its error is \(\sum x\cdot 2\,\mathrm{Re}(a^* g\, n)\big/\sum \lvert a\rvert^2 g^2\), with variance \(2N_0\int x^2 g^2 \big/ \big(\lvert a\rvert^2 (\int g^2)^2\big)\) \(= w^2/\mathrm{ENR}\). For a Gaussian PSF it is nearly efficient, because \(x\cdot g\) is exactly the derivative template. Background and box truncation bias it.
Amplitude-weighted centroid of \(\lvert y\rvert\): weights every pixel in the box with equal noise, so the error variance grows with box area (1.2–3 times the bound in the MATLAB lab).
3-point parabola or Gaussian through the peak pixel and its neighbours: cheap, but they use only three samples. As pixels shrink they see less of the PSF and lose information (about 2.2 times the bound at λ/4 and λ/8); the parabola on |y| also has a deterministic bias at coarse pixels (Module 11 showed the same for delays).

Axial: the phase is not free information. After demodulation the axial IQ PSF is the envelope times a phase \(e^{-j2kz_0}\) that also depends on \(z_0\). If that phase were known, \(\lvert\partial s/\partial z_0\rvert^2 = g'^2 + 4k^2 g^2\) (\(k = 2\pi/\lambda\)) and the axial information would be \((1 + 8k^2 w_z^2)\) times larger: 33 times (5.7 times in σ) for an axial FWHM of 0.75λ. Bubble phase depends on bubble size, shell and the exact scatterer position, so localizers treat it as unknown, and with \(a\) unknown the carrier term cancels exactly: only the envelope counts (verified in m12_nuisance_fisher.m, Part 3).

Demo 12.3

ULM localization: five localizers against the CRLB

SimulationMonte Carlo in the page (JS): Gaussian-PSF bubble images with band-limited circular noise, five localizers; bound from (12.11).

MLE (Gaussian fit, complex a)intensity centroidamplitude centroid3-point Gaussian3-point parabolaCRLB w_x/√ENR

Pixel size

Lateral FWHM1.0 λ

f₀ (for µm)

Shown ENR20 dB

Try: axial FWHM fixed at 0.75λ; noise circular complex Gaussian, white over the rectangular band bounded by the −40 dB points of each 1D PSF energy spectrum; the bubble sits anywhere within ±½ pixel of the grid centre with a random phase. At λ/4 the MLE sits on the bound (1.03–1.1×) from about 16 dB up and the intensity centroid joins it above 25 dB; the 3-point fits stay near 2× and the amplitude centroid at 1.6–2.4×. Switch to λ/2: the 3-point Gaussian improves to 1.2–1.4× (each sample now covers more of the PSF) while the parabola develops a bias floor (about 6× at 40 dB). Switch to λ/8: the bound does not move. Below about 16 dB a growing fraction of trials lock onto a noise peak (outlier fraction in the readout); one outlier in 150 trials already doubles the RMSE. Widen the PSF: everything scales with FWHM. Narrow it to 0.5λ with λ/2 pixels (pixel = FWHM): the PSF is undersampled and every localizer, the MLE included, falls off the bound.

Ultrasound link

Errico et al. (Nature, 2015) introduced ultrafast ULM and reported localization well below the diffraction limit; λ/10-type numbers in the literature are precision figures of this kind. The PALA benchmark (Heiles et al., Nature Biomedical Engineering, 2022) compared standard localizers, from weighted averages and interpolation to Gaussian fitting and radial symmetry, on simulated and in vivo data with shared metrics; read its RMSE results against the bound above, with ENR, pixel size and PSF width stated. Precision also sets how many localizations per pixel a map needs before its texture reflects vessels rather than jitter.

Do not import microscopy formulas

Single-molecule fluorescence localization (for example Thompson, Larson and Webb, Biophysical Journal, 2002) is a photon-counting problem with shot noise and pixel-integration terms. Ultrasound IQ noise is additive circular complex Gaussian, and the bubble amplitude and phase are unknown. Derive the bound from the IQ model, \(\sigma_x \ge w_x/\sqrt{\mathrm{ENR}}\), and state the noise model with every bound you quote.

12.6Nuisance parameters: what you do not know costs you

Split the parameters into the one you want, \(\theta_1\), and nuisances \(\theta_2\). The bound on \(\theta_1\) is the (1,1) entry of the inverse, which by the Schur complement is never smaller than the "everything else known" bound:

\[ \begin{aligned} \frac{1}{[\mathbf{I}^{-1}]_{11}} &= I_{11} - \mathbf{I}_{12}\,\mathbf{I}_{22}^{-1}\,\mathbf{I}_{21} \\ &\le I_{11} \\ \Rightarrow\ [\mathbf{I}^{-1}]_{11} &\ge \frac{1}{I_{11}} \end{aligned} \]

(12.12)

\([\mathbf{I}^{-1}]_{11}\): CRLB entry for the wanted parameter \(\theta_1\) with nuisances unknown
\(I_{11}\): information on the wanted parameter with every nuisance known
\(\mathbf{I}_{12},\ \mathbf{I}_{21},\ \mathbf{I}_{22}\): cross-information between wanted and nuisance parameters, and information among the nuisances

In the Gaussian case the effective information is the energy of the part of \(\partial\mathbf{s}/\partial\theta_1\) that the nuisance derivatives cannot explain: \(\lVert\mathbf{P}^{\perp}\,\partial\mathbf{s}/\partial\theta_1\rVert^2/\sigma^2\) (times 2 for complex noise). Equality holds when \(\partial\mathbf{s}/\partial\theta_1\) is orthogonal to every nuisance derivative, and the std inflation is \(1/\sqrt{1 - R^2}\).

Three ultrasound cases, all computed in m12_nuisance_fisher.m:

Unknown amplitude, width, background for a bubble centred in a window of ±3 FWHM: free (inflation 1.000), because \(g'\) is odd and \(g\), \(\partial g/\partial w\) and a constant are even. Off centre the symmetry breaks: with the bubble on the edge of a ±1.5 FWHM window (1.5 FWHM from its centre), inflation is 1.55 with amplitude unknown and 4.1 with width unknown as well (1.004 and 1.07 at 1 FWHM from the centre).
A neighbour bubble at separation \(d\) with unknown position and amplitude: inflation 3.96 at d = 0.5 FWHM, 1.38 at 1 FWHM, 1.04 at 1.5 FWHM for in-phase echoes (1.81, 1.11, 1.01 in quadrature). It does not depend on the neighbour's strength at all, only on its phase; for small \(d\) it grows as \(d^{-2}\) in phase and \(d^{-1}\) in quadrature.
Sound speed from delays: receive delays of one point target, \(\tau_e = (z + r_e)/c + t_0\), \(r_e = \sqrt{z^2 + (x_e - x_s)^2}\), 128 elements at 0.3 mm, 1 ns independent jitter (Module 11 sets that number). At z = 30 mm: \(\sigma_c\) = 0.0034 m/s if \(z\) and \(t_0\) were known, 0.065 m/s with \(z\) unknown, 0.12 m/s with \(t_0\) unknown, and 2.8 m/s with both. The intercept \(t_0\) absorbs the mean delay, depth and speed then trade along the moveout curvature, and only the small departure from a parabola separates them.

The same algebra says how to buy the information back: add data whose derivatives point in different directions (targets at several depths, several transmit angles, a calibrated lens delay \(t_0\)), or add prior information (Section 12.7).

Demo 12.4

Fisher matrix of a 1D bubble profile: the cost of each unknown

Approx. modelGaussian PSF profile (FWHM = 1) sampled at FWHM/16, circular complex white noise; Fisher matrix computed exactly from the Jacobian of this model.

|s(x)|∂s/∂x₁ (wanted)nuisance derivativespart of ∂s/∂x₁ they explain

amplitude a₁ unknown (complex) PSF width unknown background unknown (complex) neighbour bubble (x₂, a₂ unknown)

Separation d1.00 FWHM

Neighbour phase ψ0°

Neighbour amplitude |a₂/a₁|1.00

Bubble offset from window centre0.00 FWHM

Window half-width3.0 FWHM

Try: Gaussian PSF, FWHM = 1, sampled at FWHM/16, complex noise. The matrix is the Fisher matrix normalized to a correlation matrix (\(I_{ij}/\sqrt{I_{ii}I_{jj}}\)); row 1 is the wanted position. Turn the neighbour off with the bubble centred: every nuisance is free. Slide the bubble to the window edge: amplitude and width now correlate with position. Turn the neighbour on and drag d below 0.5 FWHM, then set ψ = 90°. Move the amplitude slider: the inflation does not change, because only the direction of the nuisance columns matters.

Demo 12.5

CRLB on sound speed from the delays of one point target

Approx. modelStraight-ray receive delays, 0° plane-wave transmit, independent Gaussian delay jitter per element; CRLB from the Jacobian of this delay model.

c onlyc, zc, t₀c, z, t₀your selection

∂τ/∂c∂τ/∂z∂τ/∂t₀∂τ/∂x_s

depth z unknown time offset t₀ unknown lateral position x_s unknown

Depth30 mm

Elements (pitch 0.3 mm)

Delay jitter σ_τ per element1.00 ns

Try: 0° plane-wave transmit, full receive aperture, independent Gaussian delay errors (in practice they are correlated, which is worse). The lower panel shows each Jacobian column across the aperture, scaled to unit peak: \(\partial\tau/\partial c\), \(\partial\tau/\partial z\) and \(\partial\tau/\partial t_0\) are all "constant plus a bowl", nearly collinear, which is why each added unknown costs so much. \(\partial\tau/\partial x_s\) is odd, so \(x_s\) costs nothing for a centred target. Tick \(t_0\): \(\sigma_c\) jumps by more than an order of magnitude and grows steeply with depth (F# rises, the bowl flattens).

12.7Biased estimators, shrinkage and the Bayesian view

The CRLB constrains unbiased estimators. Give up unbiasedness and the MSE can drop below \(1/I\). Shrink an unbiased estimate \(\bar y \sim \mathcal{N}(\theta, s^2)\), \(s^2\) = CRLB, toward zero:

\[ \begin{aligned} \hat\theta = c\,\bar y:\quad b &= -(1 - c)\,\theta, \quad \operatorname{var} = c^2 s^2 \\ \mathrm{MSE} &= c^2 s^2 + (1 - c)^2\theta^2 \\ \mathrm{MSE} < s^2 &\iff \lvert\theta\rvert < s\sqrt{\frac{1 + c}{1 - c}} \end{aligned} \]

(12.13)

\(\hat\theta,\ \theta\): shrinkage estimate and true parameter
\(c\): shrinkage factor, 0 to 1 (1 = unbiased)
\(\bar y\): unbiased estimate on the bound
\(b\): bias of \(\hat\theta\)
\(s^2\): variance of \(\bar y\), the CRLB
\(\mathrm{MSE},\ \operatorname{var}\): mean-squared error and variance of \(\hat\theta\)

c = 0.8: MSE = 0.64 \(s^2\) at \(\theta = 0\), break-even at \(\lvert\theta\rvert = 3s\), and worse than the bound beyond. The gain is real only where the true value is near the shrinkage target.

For \(p \ge 3\) independent means the result is stronger. The James–Stein estimator \(\hat\theta = \big(1 - (p - 2)/\lVert\mathbf{y}\rVert^2\big)\mathbf{y}\) has total MSE below \(p\) (the MLE's) for every \(\theta\), with risk exactly 2 at \(\theta = 0\) whatever \(p\) is. Pooling information across parameters is why regularized maps of many pixels can beat pixel-by-pixel estimates everywhere on average.

Bayesian estimation and MAP

Treat \(\theta\) as random with a prior \(p(\theta)\). The posterior mean minimizes the MSE averaged over the prior (MMSE), and the posterior mode is the MAP estimate. For a Gaussian prior \(\theta \sim \mathcal{N}(0, \tau^2)\) and data \(\bar y \sim \mathcal{N}(\theta, s^2)\):

\[ \begin{aligned} \mathrm{E}[\theta \mid \bar y] &= c\,\bar y, \quad c = \frac{\tau^2}{\tau^2 + s^2} \\ \mathrm{E}\big[(\hat\theta - \theta)^2\big] &= \frac{1}{1/s^2 + 1/\tau^2} \\ &= \text{Bayesian CRB} \end{aligned} \]

(12.14)

\(\mathrm{E}[\theta\mid\bar y]\): posterior mean of \(\theta\) given \(\bar y\)
\(c\): shrinkage factor, set by the prior
\(\bar y\): unbiased data estimate, \(\bar y \sim \mathcal{N}(\theta, s^2)\)
\(\tau^2\): prior variance of \(\theta\)
\(s^2\): data variance (the CRLB of \(\bar y\))
\(\hat\theta\): posterior-mean estimate
\(\text{Bayesian CRB}\): Bayesian Cramér-Rao bound (Van Trees)

Prior information adds to data information. The shrinkage factor is not a tuning knob; it is the prior. The Bayesian (Van Trees) bound is \(\mathrm{E}[(\hat\theta - \theta)^2] \ge 1/(\mathrm{E}[I(\theta)] + I_{\mathrm{prior}})\), and here the posterior mean meets it. MAP with a Gaussian prior is Tikhonov-regularized least squares (Module 13).

Demo 12.6

Shrinkage and James–Stein below the CRLB

Exact formulaSimulationTheory curves from (12.13) and the exact James–Stein risk; dots and positive-part curve from Monte Carlo in the page (JS).

MSE of c·ȳ (theory)Monte CarloCRLB (unbiased)

James–Stein risk/p (exact)positive-part JS (Monte Carlo)MLE risk/p = 1

Shrinkage factor c0.80

Number of means p10

Try: c = 1 is the unbiased estimator on the bound. Lower c: the MSE dips below the CRLB near θ = 0 and blows up for large |θ|. The readout gives the prior std τ for which this c is the posterior mean, and the Bayes MSE c·s² it earns when that prior is right. Raise p: the James–Stein curve falls further below 1 near the origin and never crosses it.

Ultrasound link

Regularized sound-speed maps, smoothed attenuation maps and learned localizers are biased estimators with a built-in prior. They can beat the CRLB in MSE when the tissue looks like the prior, and they pull the answer toward the prior when it does not (a lesion's sound speed shrunk toward the background, a bubble pulled toward the training-set PSF centre). Report bias against the true value over the range of interest, not just RMSE at one operating point.

12.8Visual summary

Each card is a question an expert answers in five seconds, the answer with its number, and the picture to keep in mind. Every curve is computed from the formula or the simulation named in the card.

12.9How experts think about it

An expert reads every precision claim as three questions. What is the data model (signal shape, noise law and covariance, real or complex, which parameters are unknown)? What does the Fisher information say for that model, roughly "energy of the signal's derivative over noise density"? And is the estimator in the regime where the local bound means anything (above threshold, unbiased, isolated targets)? Then they compare the claimed number with w/√ENR or its analogue and ask which assumption buys the difference.

MSE is bias squared plus variance, and the CRLB only bounds the variance of unbiased estimators.

\(\mathrm{MSE} = b^2 + \operatorname{var}\)biased: \(\operatorname{var} \ge (1+b')^2/I\)

Write \(\hat\theta - \theta = (\hat\theta - \mathrm{E}\hat\theta) + b\); the cross term has zero mean (12.1). An estimator with 2 µm bias and 3 µm std has RMSE \(\sqrt{13} = 3.6\) µm. Averaging more data shrinks only the 3 µm. For a biased estimator the bound becomes \((1 + b'(\theta))^2/I\), which can be far below \(1/I\) when the bias slope is negative (shrinkage has \(b' = -(1 - c)\)).

Exact formula

Says nothing about the bias itself; an estimator can have tiny variance and a large bias (pixel-centre localization, clipped search ranges).

Fisher information is the curvature of the log-likelihood: a sharp peak means a small error.

\(I = -\mathrm{E}[\partial^2 \ln p]\)\(\sigma \ge 1/\sqrt{I}\)independent data add

Near its maximum the log-likelihood is approximately \(-\tfrac12 I\,(\theta - \hat\theta)^2\), so its width is \(1/\sqrt{I}\). For the mean of N samples of variance \(\sigma^2\), \(I = N/\sigma^2\): four times the data, four times the curvature, half the standard deviation. Information from independent measurements (frames, elements, angles) simply adds.

Exact formula

Expected curvature; one realization can be flatter. Multimodal likelihoods (RF correlation, periodograms at low SNR) have several peaks, and the curvature of one says nothing about jumping to another.

In Gaussian noise, Fisher information is the energy of the signal's derivative divided by the noise level.

real: \((1/\sigma^2)\sum(\partial s/\partial\theta)^2\)complex: \((2/\sigma^2)\,\mathrm{Re}\sum\lvert\partial s/\partial\theta\rvert^2\)

For \(\mathbf{s}(\theta) + \mathbf{w}\) with \(\mathbf{w} \sim \mathcal{N}(0, \sigma^2\mathbf{I})\), \(\ln p = -\lVert\mathbf{y} - \mathbf{s}\rVert^2/(2\sigma^2) + \text{const}\); differentiating twice and taking the expectation removes the residual term and leaves \(\mathbf{J}^{T}\mathbf{J}/\sigma^2\) (12.6). Everything else in this module is this formula with a different \(\mathbf{J}\): PSF derivative for position, \(n\cdot s\) for frequency, \(\dot s\) for delay, moveout derivatives for sound speed.

Exact formula

Deterministic signal in additive Gaussian noise. For random signals (speckle, blood), the parameter also changes the signal covariance and a different (Whittle or covariance-based) formula applies.

Complex IQ noise carries a factor 2: count \(\sigma^2\) per quadrature, not per sample.

\(\mathrm{E}\lvert w\rvert^2 = \sigma^2 \to \sigma^2/2\) per quadrature\(\operatorname{var}(\hat a) \ge \sigma^2/(2\lVert s\rVert^2)\)

A complex sample is two real samples, each with variance \(\sigma^2/2\), so the real formula gives \((2/\sigma^2)\sum\lvert\partial s\rvert^2\) (12.7). A phase-known amplitude from IQ thus has half the variance of the same pulse in real noise of variance \(\sigma^2\) (the MATLAB lab gives ratios 0.99, 1.00, 0.99 against \(\sigma^2/(2\lVert s\rVert^2)\), \(\sigma^2/\lVert s\rVert^2\), \(\sigma^2/\lVert s\rVert^2\)). The same factor separates Kay's real-sinusoid bound from the complex-tone bound at equal SNR.

Simulation300 circular complex noise samples (JS, fixed seed) and the Gaussian pdf per quadrature (std σ/√2) against real noise of std σ.

Requires circular noise (equal, uncorrelated quadratures). Demodulation with a poor low-pass or I/Q imbalance breaks circularity; then use the 2N-dimensional real model.

Localization precision is PSF width over the square root of SNR: \(\sigma_x = w_x/\sqrt{\mathrm{ENR}}\).

w = FWHM/2.35520 dB → FWHM/23.530 dB → FWHM/74.5

For a Gaussian PSF, \(\int g'^2/\int g^2 = 1/(2w^2)\), so \(I = (2\lvert a\rvert^2/N_0)\cdot E_g/(2w^2) = \mathrm{ENR}/w^2\) (12.11). With FWHM = λ = 308 µm (5 MHz) and ENR = 30 dB, \(\sigma_x\) = 4.1 µm. λ/10 needs only 12.6 dB of ENR in principle, which is why sub-λ/10 precision claims are physically easy and say little about the image.

Approx. modelGaussian PSF; pdfs drawn at the CRLB.

Isolated bubble, Gaussian PSF, noise white over the signal band, amplitude and phase unknown, above threshold. Real PSFs have sidelobes and depth-dependent shape; motion and tracking add error.

λ/10 in ULM is localization precision, not resolution of two simultaneous bubbles.

×3.96 at d = 0.5 FWHM×1.38 at 1 FWHM≈1 beyond 1.5 FWHM

Add a neighbour with unknown position and amplitude: the derivative of bubble 1 overlaps the nuisance columns of bubble 2, and the std inflation is \(1/\sqrt{1 - R^2}\) (12.12). For in-phase echoes it is 3.96 at 0.5 FWHM separation and grows as \(d^{-2}\) (quadrature: 1.81, growing as \(d^{-1}\)). That is why ULM needs sparse bubbles, and why high-concentration methods become sparse-recovery problems.

Approx. modelGaussian PSF; inflation from the numerical Fisher matrix.

Fisher (local) analysis assumes the number of bubbles is known; deciding whether there are one or two is a detection problem that fails earlier.

Pixels finer than the noise correlation length add no information; they only help the estimator.

bound depends on \(N_0\) and PSFiid-per-pixel sim: \(\sigma \propto p\) (fake)

Beamformed noise is band-limited: oversampled pixels carry correlated noise, and the continuous bound \(w/\sqrt{\mathrm{ENR}}\) is reached by a matched-filter MLE at λ/2, λ/4 or λ/8 alike (MATLAB: 0.95–1.03 times the bound). A simulation adding independent noise per pixel at fixed per-pixel SNR makes \(\mathrm{ENR} \propto 1/p^2\), inventing a factor-2 gain in σ per halving of the pixel.

Approx. modelScaling law: iid per-pixel noise at fixed pixel SNR gives ENR ∝ 1/p².

Pixels coarser than the Nyquist spacing of the PSF band do lose information (aliasing), and every simple estimator is pixel-dependent.

For a Gaussian PSF the intensity-weighted centroid is nearly efficient; amplitude centroids and 3-point fits are not.

IC: \(\operatorname{var} \to w^2/\mathrm{ENR}\)3-point fits ≈ 2.2× at λ/4

At high ENR, \(\lvert y\rvert^2 \approx \lvert a\rvert^2 g^2 + 2\,\mathrm{Re}(a^* g\, n)\), so the intensity centroid's error is \(\sum x\cdot 2\,\mathrm{Re}(a^* g\, n)/\sum\lvert a\rvert^2 g^2\), with variance \(2N_0\int x^2 g^2/(\lvert a\rvert^2(\int g^2)^2)\) \(= w^2/\mathrm{ENR}\) because \(x\cdot g \propto g'\) is the optimal template. Amplitude weights add noise from the whole box; 3-point fits discard most of the PSF.

Approx. modelGaussian PSF, high-ENR linearization.

Background, box truncation and neighbour bubbles bias the centroid; at low ENR the \(\lvert n\rvert^2\) term adds error (1.3–1.5× at 16–20 dB in the lab).

The CRLB is a local bound: check the threshold in integrated SNR before quoting it.

N = 32 tone: on bound ≥ 13 dB \(N\cdot\mathrm{SNR}\)≈180× at 11 dB

The bound comes from the curvature at the true peak. Once a noise peak somewhere else can win, errors become a mixture of small jitter and large jumps, and the MSE heads to the no-information level (\(\pi^2/3\) for frequency). The knee depends on the energy in the whole record, \(N\cdot\mathrm{SNR}\) or ENR, not the per-sample SNR. Ziv–Zakai and Barankin bounds include those jumps and reproduce the knee.

SimulationMonte Carlo, values from m12_frequency_crlb.m (Octave), N = 32.

The exact threshold depends on the search range (number of competing cells) and on the estimator; measure it by simulation for your geometry.

Frequency precision improves as \(N^{-3/2}\): the ensemble length is a lever arm.

\(\operatorname{var}(\hat\omega) \ge 6/(\mathrm{SNR}\cdot N(N^2-1))\)\(\sigma_v = c\cdot\mathrm{PRF}\cdot\sigma_\omega/(4\pi f_0)\)

A frequency error shows up as a phase ramp, so samples far from the record centre carry more information (\(I \propto \sum(n - \bar n)^2\), (12.8)). Doubling N from 16 to 32 cuts \(\sigma_\omega\) by 2.83. M = 16, 10 dB, PRF 5 kHz, \(f_0\) 5 MHz: \(\sigma_\omega\) = 0.0121 rad/sample, 9.65 Hz, 1.49 mm/s.

Exact formula

Deterministic constant tone. Blood is a random process with spectral width and acceleration; longer ensembles also see more tissue motion and decorrelation.

Lag-one autocorrelation wastes information at high SNR: \(N(N+1)/(6(N-1))\) times the bound in variance.

5.7× at N = 32\(\operatorname{var} \to 1/(\mathrm{SNR}\,(N-1)^2)\)

At high SNR the phase of \(\sum y[n+1]\,y^*[n]\) is the average of \(N - 1\) phase differences, which telescopes to \((\varphi_{N-1} - \varphi_0)/(N - 1)\): only the two end samples count. Variance \(2\cdot(1/(2\,\mathrm{SNR}))/(N - 1)^2\) \(= 1/(\mathrm{SNR}\,(N - 1)^2)\), against \(6/(\mathrm{SNR}\cdot N(N^2 - 1))\). At moderate SNR it is worse (about 15× at 10 dB, N = 32).

Exact formula

For broad Doppler spectra (real blood) the tone MLE is not the right reference, and lag-one remains a robust, cheap choice.

A nuisance parameter never helps: \([\mathbf{I}^{-1}]_{ii} \ge 1/I_{ii}\), with equality when its derivative is orthogonal.

inflation \(= 1/\sqrt{1 - R^2}\)symmetry decouples

The Schur complement removes from \(\partial s/\partial\theta_1\) whatever the nuisance derivatives can mimic. For a centred symmetric PSF, \(g'\) is odd while \(g\), \(\partial g/\partial w\) and a constant are even, so unknown amplitude, width and background cost nothing (inflation 1.000). At the edge of a ±1.5 FWHM window the same unknowns cost 1.55× (amplitude) and 4.1× (with width).

Approx. modelGaussian PSF.

Local analysis; nuisances that make the problem nonlinear or multimodal (unknown number of targets) cost more than the Fisher matrix shows.

With unknown phase, only the envelope localizes axially; a known carrier phase would be worth \(1 + 8k^2 w_z^2\) in information.

33× info, 5.7× σ (\(\mathrm{FWHM}_z\) = 0.75λ)unknown \(a\) ⇒ carrier cancels

The demodulated axial PSF is \(g(z - z_0)\,e^{-j2kz_0}\), \(k = 2\pi/\lambda\) (the 2 is the round trip). If the phase is tied to \(z_0\), \(\partial s/\partial z_0 = -(g' + j2kg)\,e^{-j2kz_0}\), so \(\lvert\partial s/\partial z_0\rvert^2 = g'^2 + 4k^2 g^2\) integrates to \((1/(2w^2) + 4k^2)\cdot E\): a ratio \(1 + 8k^2 w^2 = 33\) for \(w_z = 0.75\lambda/2.355 = 0.32\lambda\). With a complex amplitude unknown, the phase column absorbs the carrier term exactly (Fisher ratio 1.0000 in the lab). This is the localization version of RF-versus-envelope delay estimation in Module 11.

Approx. modelGaussian axial envelope, FWHM 0.75λ.

Phase tracking across frames (for a bubble that moves slowly) can recover part of that information, but phase wraps every λ/4 of axial motion.

Sound speed from one target is poor once the time offset and depth are unknown: they trade against c.

z = 30 mm, 1 ns: 0.065 m/s (\(c, z\))2.8 m/s (\(c, z, t_0\))

\(\tau_e = (z + r_e)/c + t_0\). In the small-aperture limit the moveout is a parabola: intercept \(2z/c\) (merged with \(t_0\)) and curvature \(1/(2zc)\). With \(t_0\) known, intercept and curvature give \(c\) (\(\sigma_c \approx \sqrt{180}\,\sigma_\tau z c^2/(\sqrt{N} D^2)\) = 0.057 m/s at 30 mm against 0.065 exact). With \(t_0\) unknown only the departure from a parabola separates \(c\) and \(z\), and \(\sigma_c\) rises 40-fold and grows roughly as \(z^3\) at depth.

Approx. modelStraight-ray receive delays, 128 elements at 0.3 mm, z = 30 mm.

Independent delay errors and a known geometry; correlated (aberration) delay errors and refraction make it worse. Multiple targets, angles and a calibrated lens delay break the degeneracy.

Least squares is the MLE in white Gaussian noise; with coloured noise, whiten first.

\(\hat\theta = \arg\min\,(\mathbf{y}-\mathbf{s})^{H}\mathbf{C}^{-1}(\mathbf{y}-\mathbf{s})\)linear model: LS efficient at any N

The Gaussian log-likelihood is a quadratic form in the residual, so ML is (weighted) least squares (12.2). In the linear model \(\operatorname{cov}(\hat\theta) = \sigma^2(\mathbf{H}^{T}\mathbf{H})^{-1}\) equals the CRLB exactly. Ignoring correlation (ordinary LS on beamformed pixels, adjacent-element delays, overlapping spectral gates) leaves the estimate unbiased but inefficient, and makes its reported error bars wrong.

Simulation260 correlated noise pairs (ρ = 0.9, JS, fixed seed) and the same samples after whitening by \(\mathbf{C}^{-1/2}\).

Heavy-tailed noise or outliers (clutter spikes, false peaks): Gaussian ML is not robust; use a robust loss or reject outliers first.

The MLE is invariant, not unbiased: \(\lvert\hat a\rvert^2\) overestimates the power by the noise term.

\(\mathrm{E}\lvert\hat a\rvert^2 = \lvert a\rvert^2 + \sigma^2/\lVert s\rVert^2\)bias \(= 1/\mathrm{ENR}\) relative

Invariance makes \(\lvert\hat a\rvert^2\) the MLE of \(\lvert a\rvert^2\), but \(\mathrm{E}\lvert\hat a\rvert^2 = \lvert a\rvert^2 + \mathrm{E}\lvert\hat a - a\rvert^2\) adds the estimator variance (0.0474 measured vs 0.0470 predicted in the lab). At 10 dB ENR the power is 10 % high. The same mechanism biases Power Doppler, coherence and envelope statistics at low SNR; subtract a noise-power estimate to correct it.

Exact formula

Asymptotically (high ENR) the bias vanishes, as ML theory promises.

Average intensity, not decibels: the dB average is biased by −2.51 dB and is 61 % efficient.

\(\mathrm{E}[\ln I] = \ln\mu - \gamma\)\(\operatorname{var}\ln I = \pi^2/6 = 1.64\)

Speckle intensity is exponential with mean \(\mu\); its sample mean is the efficient MLE (\(\operatorname{var} = \mu^2/N\), i.e. \(1/N\) for \(\ln\mu\)). The log of an exponential variable has mean \(\ln\mu - 0.577\) (−2.51 dB in \(10\log_{10}\)) and variance \(\pi^2/6\), so averaging dB values estimates \(\ln\mu\) with \(1.64/N\): efficiency \(6/\pi^2 = 0.61\).

Exact formula

Fully developed speckle (exponential intensity). With few scatterers or a coherent component the bias changes, which is itself an estimator of the envelope statistics (Module 26).

A biased estimator can beat the CRLB in MSE, but only where its prior is right.

\(c\bar y\) break-even \(\lvert\theta\rvert = s\sqrt{(1+c)/(1-c)}\)JS risk = 2 at \(\theta = 0\)

Shrinkage trades a bias \((1 - c)\theta\) for variance \(c^2 s^2\): c = 0.8 gives 0.64 of the CRLB at \(\theta = 0\) and breaks even at \(3s\) (12.13). It is exactly the posterior mean for a prior std \(\tau = s\sqrt{c/(1 - c)}\), whose Bayes MSE \(1/(1/s^2 + 1/\tau^2) = c\,s^2\) is the Bayesian CRB (12.14). James–Stein (\(p \ge 3\)) improves total MSE for every \(\theta\).

Exact formula

A wrong prior makes the estimator worse than unbiased far from the target (MSE 2.4× the CRLB at \(\theta = 3\tau\) for \(\tau = 1.5s\)); learned estimators inherit their training distribution as the prior.

12.10Cheat sheet

Concept	Formula	Ultrasound use
MSE decomposition	\(\mathrm{MSE} = b^2 + \operatorname{var}\)	Report bias and jitter of any localizer or tracker separately
Fisher information	\(I = \mathrm{E}[(\partial\ln p/\partial\theta)^2] = -\mathrm{E}[\partial^2\ln p/\partial\theta^2]\)	Precision limit of any parameter
CRLB (scalar, vector)	\(\operatorname{var}\hat\theta \ge 1/I\); \(\operatorname{var}\hat\theta_i \ge [\mathbf{I}^{-1}]_{ii}\); \(\operatorname{var}\hat\alpha \ge \nabla f^{T}\mathbf{I}^{-1}\nabla f\)	Error bars on SoS, attenuation, delays, positions
Gaussian, real noise	\(\mathbf{I} = \mathbf{J}^{T}\mathbf{C}^{-1}\mathbf{J}\) (white: \(\mathbf{J}^{T}\mathbf{J}/\sigma^2\))	RF data
Gaussian, complex noise	\(\mathbf{I} = 2\,\mathrm{Re}\{\mathbf{J}^{H}\mathbf{C}^{-1}\mathbf{J}\}\) (white, \(\mathrm{E}\lvert w\rvert^2 = \sigma^2\): \(2\,\mathrm{Re}\,\mathbf{J}^{H}\mathbf{J}/\sigma^2\))	IQ data
Amplitude of known pulse (IQ)	real \(a\): \(\sigma^2/(2\lVert s\rVert^2)\); complex \(a\): \(\mathrm{E}\lvert\hat a - a\rvert^2 \ge \sigma^2/\lVert s\rVert^2\)	Echo amplitude, BSC
Concentrated ML, unknown complex \(a\)	\(\hat\theta = \arg\max\,\lvert\mathbf{g}(\theta)^{H}\mathbf{y}\rvert^2/\lVert\mathbf{g}(\theta)\rVert^2\)	Periodogram, matched-filter localization, correlation delay
Complex tone frequency	\(\operatorname{var}\hat\omega \ge 6/(\mathrm{SNR}\cdot N(N^2-1))\)	Slow-time Doppler; \(\sigma_v = c\cdot\mathrm{PRF}\cdot\sigma_\omega/(4\pi f_0)\)
Real sinusoid (Kay Ex. 3.14)	\(\operatorname{var}\hat f \ge 12/((2\pi)^2\eta N(N^2-1))\), \(\eta = A^2/(2\sigma^2)\)	Real-valued tone fits
Lag-one (Kasai) efficiency	high SNR: \(N(N+1)/(6(N-1))\) × CRLB	5.7× at N = 32
ULM localization, Gaussian PSF	\(\sigma_x \ge w_x/\sqrt{\mathrm{ENR}}\), \(w = \mathrm{FWHM}/2.355\), \(\mathrm{ENR} = \lVert s\rVert^2/N_0\)	Precision of bubble positions; real noise: \(\sqrt{2}\,w/\sqrt{\mathrm{ENR}_r}\)
Peak-pixel SNR vs ENR	\(\mathrm{SNR}_{\mathrm{pk}} \approx \mathrm{ENR} - 4.7\ldots5\ \mathrm{dB}\) (noise white over the rectangle bounded by the −40 dB points of each 1D PSF spectrum; 3.6 dB for the elliptical −40 dB band)	Convert paper SNRs to ENR
Known vs unknown carrier phase	\(I_{\mathrm{known}}/I_{\mathrm{env}} = 1 + 8k^2 w_z^2\), \(k = 2\pi/\lambda\)	33× for \(\mathrm{FWHM}_z\) = 0.75λ
Nuisance cost	\(1/[\mathbf{I}^{-1}]_{11} = I_{11} - \mathbf{I}_{12}\mathbf{I}_{22}^{-1}\mathbf{I}_{21}\); inflation \(1/\sqrt{1-R^2}\)	Neighbour bubbles, \(t_0\) and \(z\) in SoS fits
Shrinkage	\(\mathrm{MSE} = c^2 s^2 + (1-c)^2\theta^2\)	Regularized maps; break-even \(s\sqrt{(1+c)/(1-c)}\)
Bayesian CRB	\(\mathrm{E}(\hat\theta-\theta)^2 \ge 1/(\mathrm{E}[I] + I_{\mathrm{prior}})\)	MAP / Tikhonov, learned estimators
Intensity averaging	\(\mathrm{E}[10\log_{10} I] = 10\log_{10}\mu - 2.51\ \mathrm{dB}\)	Average intensity, not dB

12.11Interview and qualifying-exam questions

I1What is the Cramér–Rao bound, and what does it assume? · screen

A lower bound on the variance of any unbiased estimator: \(\operatorname{var}(\hat\theta) \ge 1/I(\theta)\), \(I = \mathrm{E}[(\partial\ln p/\partial\theta)^2]\). It assumes the support of \(p(y;\theta)\) does not depend on \(\theta\) and the log-likelihood is differentiable; it is local (curvature at the true value) and says nothing about bias. It is reached exactly in the linear Gaussian model and asymptotically by the MLE. Follow-up: "And for biased estimators?" \((1 + b')^2/I\), and the MSE can be lower than \(1/I\).

I2A paper claims λ/20 resolution with ULM. What do you ask? · screen

Whether λ/20 is localization precision (std of an isolated bubble's position) or the separability of two bubbles. Then the ENR or SNR definition and noise model (\(\sigma_x = \mathrm{FWHM}/(2.355\sqrt{\mathrm{ENR}})\): λ/20 with FWHM ≈ λ needs only ENR ≈ 18.6 dB), whether outliers were excluded, how motion and tracking errors were handled, and how resolution was measured on the final map (for example a Fourier-ring-correlation split of the localizations, or the smallest separated vessel pair). Follow-up: "Why is precision so cheap?" Because the bound scales as \(1/\sqrt{\mathrm{ENR}}\) and bubbles are bright.

I3Derive the localization CRLB for a Gaussian PSF in complex white noise. · onsite

\(s = a\,g(x - x_0)\,g(z - z_0)\), noise density \(N_0\). \(I_{x_0x_0} = (2/N_0)\lvert a\rvert^2\iint\lvert\partial g/\partial x\rvert^2\) \(= (2\lvert a\rvert^2/N_0)\cdot\pi w_x w_z/(2w_x^2)\) \(= \mathrm{ENR}/w_x^2\). Cross terms with \(\mathrm{Re}\,a\), \(\mathrm{Im}\,a\) and \(z_0\) vanish (odd × even), so \(\sigma_x \ge w_x/\sqrt{\mathrm{ENR}}\). Example: FWHM 1λ at 5 MHz, ENR 30 dB: 308/(2.355·31.6) = 4.1 µm. Follow-up: "Real RF instead?" \(\sqrt{2}\,w/\sqrt{\mathrm{ENR}_r}\) with \(\mathrm{ENR}_r = \sum s^2/\sigma^2\), and the axial bound changes if the carrier phase is known.

I4Where does the factor 2 between real and complex noise come from? · onsite

With \(\mathrm{E}\lvert w\rvert^2 = \sigma^2\), each quadrature has \(\sigma^2/2\); summing the real-data information of both quadratures gives \(I = (2/\sigma^2)\,\mathrm{Re}\sum\lvert\partial s\rvert^2\). Phase-known amplitude from IQ: \(\sigma^2/(2\lVert s\rVert^2)\); real pulse in real noise \(\sigma^2\): \(\sigma^2/\lVert s\rVert^2\). Same for frequency: complex tone \(6/(\mathrm{SNR}\cdot N(N^2-1))\) vs Kay's real \(12/(\eta N(N^2-1))\). Follow-up: "Does demodulating RF to IQ change the bound?" No; it is an invertible map of band-limited data. Only the bookkeeping of \(\sigma^2\) changes.

I5Why can an intensity-weighted centroid be nearly efficient while an amplitude-weighted centroid is not? · onsite

Linearize at high ENR: \(\lvert y\rvert^2 \approx \lvert a\rvert^2 g^2 + 2\,\mathrm{Re}(a^* g\, n)\), so the intensity centroid's error is a noise projection onto \(x\cdot g\), which is proportional to \(g'\), the optimal (matched) template: variance \(w^2/\mathrm{ENR}\), the bound. With \(\lvert y\rvert\) weights the error is a projection onto \(x\cdot 1\) over the box: equal-weight noise from every pixel, so variance grows with box area. Follow-up: "What breaks the intensity centroid?" Background (pulls toward the box centre), truncation, neighbours, and the \(\lvert n\rvert^2\) term at low ENR.

I6How does Doppler frequency precision scale with ensemble length, and where is the threshold? · onsite

\(\operatorname{var}(\hat\omega) \ge 6/(\mathrm{SNR}\cdot N(N^2 - 1))\), so \(\sigma \propto N^{-3/2}\); N = 16 → 32 cuts σ by 2.83. The periodogram MLE follows it down to an integrated \(N\cdot\mathrm{SNR}\) of about 13 dB (N = 32); 2 dB lower it is two orders of magnitude worse because noise bins win. Lag-one is 5.7× in variance above the bound at N = 32 even at high SNR. Follow-up: "Convert to velocity": \(\sigma_v = c\cdot\mathrm{PRF}\cdot\sigma_\omega/(4\pi f_0)\): 1.49 mm/s for N = 16, 10 dB, PRF 5 kHz, 5 MHz.

I7How would you bound sound-speed estimation from delay data, and what dominates the bound? · research

Write the delay model \(\tau_e(c, z, x_s, t_0)\), its Jacobian, \(\mathbf{I} = \mathbf{J}^{T}\mathbf{C}_\tau^{-1}\mathbf{J}\) with the delay covariance from the Walker–Trahey/Carter bound (Module 11), and take \([\mathbf{I}^{-1}]_{cc}\). The nuisances dominate: for one target at 30 mm, 128 elements, 1 ns jitter, \(\sigma_c\) goes from 0.065 m/s (\(z\) unknown) to 2.8 m/s (\(z\) and \(t_0\) unknown), because \(t_0\), \(z\) and \(c\) all produce nearly the same bowl-shaped moveout. Break it with several depths and angles, a calibrated lens delay, or a prior. Follow-up: "Correlated delay errors?" Use the full \(\mathbf{C}_\tau\); aberration-like errors are low-order and project onto the same bowl, so they hurt more than white jitter.

I8A learned localizer beats the CRLB on your test set. Is that suspicious? · research

Not necessarily: the CRLB bounds unbiased estimators. A network trained on a distribution of PSFs and positions acts as a Bayesian estimator with that prior, and the relevant bound is the Bayesian CRB \(1/(\mathrm{E}[I] + I_{\mathrm{prior}})\). Check the bias map against true sub-pixel position (pull toward pixel centres or training positions), test out of distribution (other PSF widths, depths, neighbours), and compare the MSE against the Bayesian bound for the same prior. Follow-up: "Can any estimator beat the Bayesian bound?" Not in Bayes MSE under regularity conditions.

I9Why is the MLE of echo power from IQ biased, and how do you correct it? · onsite

Invariance makes \(\lvert\hat a\rvert^2\) the MLE of \(\lvert a\rvert^2\), but \(\mathrm{E}\lvert\hat a\rvert^2 = \lvert a\rvert^2 + \sigma^2/\lVert s\rVert^2\): the estimator variance adds. At ENR = 10 dB the power is 10 % (0.41 dB) high. Subtract an estimate of \(\sigma^2/\lVert s\rVert^2\) from a noise-only region or noise frames. Follow-up: "Where else?" Power Doppler, coherence factor and correlation coefficients at low SNR, all biased by the noise term.

I10Two bubbles are 0.5 FWHM apart. What happens to the precision, and what does it mean for high-concentration ULM? · research

With both positions and amplitudes unknown, the Fisher matrix couples them: std inflation 3.96 for in-phase echoes (1.81 in quadrature), rising as \(d^{-2}\) (\(d^{-1}\)) as they approach. Before that, a single-bubble localizer returns a biased "mid-point". High concentration therefore needs model-order selection (one bubble or two), sparse recovery (Module 13), or temporal separation (separable bubble dynamics). Follow-up: "Does amplitude ratio matter?" Not in the Fisher inflation (only the relative phase), but a weak neighbour is harder to detect, which fails first.

I11To estimate mean backscattered power from speckle, do you average intensity or dB values? · screen

Intensity. The sample mean of exponential intensities is the efficient MLE. Averaging dB values gives a biased result (−2.51 dB, Euler's \(\gamma\) in \(10\log_{10}\)) with variance \(\pi^2/6 = 1.64\) times the bound for \(\ln\mu\): 61 % efficiency, i.e. the same precision needs 64 % more independent samples. Follow-up: "How many independent samples in a gate?" Count speckle cells, not pixels (Module 07).

12.12MATLAB lab

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

Download the scripts: m12_ulm_localization_crlb.m · m12_crlb_amplitude_mle.m · m12_frequency_crlb.m · m12_nuisance_fisher.m · m12_shrinkage_mse.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m12_ulm_localization_crlb.mMLE, intensity and amplitude centroids, 3-point Gaussian and parabola vs the CRLB over ENR and pixel size (λ/2, λ/4, λ/8) · change fwhmX, pixList, enrList · runtime ≈ 15 s

%% m12_ulm_localization_crlb.m
% Purpose : ULM localization precision for one isolated bubble in a 2D
%           beamformed IQ image, compared with the Cramer-Rao bound.
%   Model : s(x,z) = a * gx(x - x0) * gz(z - z0), Gaussian PSF with
%           lateral/axial FWHM, a complex and unknown (random phase).
%           Noise: circular complex Gaussian, white inside the system band
%           |kx| <= 0.483/wx, |kz| <= 0.483/wz (PSF energy spectrum -40 dB),
%           zero outside, density N0. This mimics beamformed noise, which
%           is band-limited: pixels finer than Nyquist add no information.
%   Bound : sigma_x >= wx / sqrt(ENR), sigma_z >= wz / sqrt(ENR),
%           w = FWHM / 2.3548, ENR = ||s||^2 / N0 = |a|^2 pi wx wz / N0
%           (integrated SNR; complex noise; amplitude and phase unknown).
%   Estimators (all from the pixel with the largest |IQ|):
%           WA  : amplitude-weighted centroid of |IQ| in a +-FWHM box
%           IC  : intensity-weighted centroid of |IQ|^2 in the same box
%           PAR : separable 3-point parabola on |IQ|
%           GAU : separable 3-point Gaussian fit (parabola on log|IQ|)
%           MLE : maximize |<y, g(x,z)>|^2 / ||g||^2 (Newton, 2D)
% Change  : fwhmX, fwhmZ, pixList, enrList, nTrial, f0.
% Expect  : MLE within 0.95-1.05x of the bound for ENR >= 20 dB at every
%           pixel size; IC (intensity centroid) 1.0-1.3x; GAU 1.2-1.4x at
%           lambda/2 but ~2.2x at lambda/4 and lambda/8 (a 3-point fit uses
%           less of the PSF as pixels shrink); PAR has a bias floor at
%           lambda/2 (5.7x at 40 dB); WA 1.2-3x; NS (pixel centre) floors at
%           pix/sqrt(12). Peak-pixel SNR = ENR - 4.8 dB (lambda/4, lambda/8).
%           At 16 dB 1% outliers (a noise peak beats the bubble) already give
%           2-2.5x; at 10 dB about half the trials are outliers (out%).
rng(1);
% 3-point vertex offset in pixels, clamped to +-0.5 (exact bound when the
% centre sample is the maximum; guards against edge peaks at low SNR)
off3 = @(v) max(-0.5, min(0.5, 0.5*(v(1) - v(3))/min(v(1) - 2*v(2) + v(3), -eps)));

%% Parameters (lengths in wavelengths)
f0 = 5e6; c = 1540; lam = c/f0;            % 5 MHz: lambda = 308 um
fwhmX = 1.0; fwhmZ = 0.75;                 % PSF FWHM (lambda), lateral and axial
wx = fwhmX/(2*sqrt(2*log(2))); wz = fwhmZ/(2*sqrt(2*log(2)));
pixList = [1/2 1/4 1/8];                   % square pixel size (lambda)
enrList = [10 16 20 25 30 40];            % integrated SNR (dB)
nTrial = 300;
names = {'MLE', 'IC ', 'WA ', 'GAU', 'PAR', 'NS '};
nEst = numel(names);
crlbX = wx ./ sqrt(10.^(enrList/10));
crlbZ = wz ./ sqrt(10.^(enrList/10));
rmseX = zeros(numel(pixList), numel(enrList), nEst);
rmseZmle = zeros(numel(pixList), numel(enrList));
outMLE = zeros(numel(pixList), numel(enrList));   % fraction of MLE outliers
snrPkOff = zeros(1, numel(pixList));             % ENR - peak-pixel SNR (dB)
tic;
for ip = 1:numel(pixList)
    pix = pixList(ip);
    nx = ceil(3*fwhmX/pix); nz = ceil(3*fwhmZ/pix);
    xs = (-nx:nx)*pix; zs = (-nz:nz)'*pix;            % pixel grid
    Nx = numel(xs); Nz = numel(zs);
    Kx = ceil(fwhmX/pix); Kz = ceil(fwhmZ/pix);       % centroid half-box (pixels)
    Lx = 2^nextpow2(2*Nx); Lz = 2^nextpow2(2*Nz);
    kx = [0:Lx/2, -Lx/2+1:-1]/(Lx*pix); kz = [0:Lz/2, -Lz/2+1:-1]'/(Lz*pix);
    mask = (abs(kz) <= 0.483/wz) * (abs(kx) <= 0.483/wx);   % rectangular band
    Eg = pi*wx*wz;                                    % integral of |g|^2 (a = 1)
    se = zeros(numel(enrList), nEst); seZ = zeros(numel(enrList), 1); nOut = zeros(numel(enrList), 1);
    sig2pix = mean(mask(:))/pix^2;                    % per-pixel noise variance (N0 = 1)
    snrPkOff(ip) = 10*log10(Eg*sig2pix);              % ENR = peak SNR * Eg * sig2pix
    for t = 1:nTrial
        x0 = (rand - 0.5)*pix; z0 = (rand - 0.5)*pix; ph = 2*pi*rand;
        Wn = (randn(Lz, Lx) + 1j*randn(Lz, Lx)) * sqrt(1/(2*pix*pix));   % density N0 = 1
        Wn = ifft2(fft2(Wn) .* mask);
        Wn = Wn(1:Nz, 1:Nx);
        G = exp(-(zs - z0).^2/(2*wz^2)) * exp(-(xs - x0).^2/(2*wx^2));
        for ie = 1:numel(enrList)
            A = sqrt(10^(enrList(ie)/10)/Eg);
            Y = A*exp(1j*ph)*G + Wn;
            M = abs(Y);
            [~, imax] = max(M(:)); [iz, ix] = ind2sub(size(M), imax);
            iz = min(max(iz, 2), Nz - 1); ix = min(max(ix, 2), Nx - 1);
            rz = max(iz - Kz, 1):min(iz + Kz, Nz); rx = max(ix - Kx, 1):min(ix + Kx, Nx);
            B = M(rz, rx); XX = ones(numel(rz), 1)*xs(rx); ZZ = zs(rz)*ones(1, numel(rx));
            wa = [sum(B(:).*XX(:)), sum(B(:).*ZZ(:))]/sum(B(:));
            ic = [sum(B(:).^2.*XX(:)), sum(B(:).^2.*ZZ(:))]/sum(B(:).^2);
            mx = M(iz, ix-1:ix+1); mz = M(iz-1:iz+1, ix).';
            par = [xs(ix) + pix*off3(mx), zs(iz) + pix*off3(mz)];
            gau = [xs(ix) + pix*off3(log(mx)), zs(iz) + pix*off3(log(mz))];
            % MLE: Newton on L(x,z) = ln|gz' Y gx|^2 - ln(||gx||^2 ||gz||^2), start at IC
            p = ic; p0 = p;
            for it = 1:8
                dx = xs.' - p(1); dz = zs - p(2);
                gx = exp(-dx.^2/(2*wx^2)); gx1 = gx.*dx/wx^2; gx2 = gx.*(dx.^2/wx^4 - 1/wx^2);
                gz = exp(-dz.^2/(2*wz^2)); gz1 = gz.*dz/wz^2; gz2 = gz.*(dz.^2/wz^4 - 1/wz^2);
                YX = Y*[gx gx1 gx2];                       % Nz x 3
                Z = gz.'*YX(:,1); Zx = gz.'*YX(:,2); Zxx = gz.'*YX(:,3);
                Zz = gz1.'*YX(:,1); Zzz = gz2.'*YX(:,1); Zxz = gz1.'*YX(:,2);
                P = abs(Z)^2;
                Px = 2*real(conj(Z)*Zx); Pz = 2*real(conj(Z)*Zz);
                Pxx = 2*(abs(Zx)^2 + real(conj(Z)*Zxx)); Pzz = 2*(abs(Zz)^2 + real(conj(Z)*Zzz));
                Pxz = 2*real(conj(Zx)*Zz + conj(Z)*Zxz);
                ex = sum(gx.^2); ex1 = 2*sum(gx.*gx1); ex2 = 2*sum(gx1.^2 + gx.*gx2);
                ez = sum(gz.^2); ez1 = 2*sum(gz.*gz1); ez2 = 2*sum(gz1.^2 + gz.*gz2);
                grad = [Px/P - ex1/ex; Pz/P - ez1/ez];
                H = [Pxx/P - (Px/P)^2 - ex2/ex + (ex1/ex)^2, Pxz/P - Px*Pz/P^2; ...
                     Pxz/P - Px*Pz/P^2, Pzz/P - (Pz/P)^2 - ez2/ez + (ez1/ez)^2];
                if H(1,1) >= 0 || det(H) <= 0, break; end
                st = -(H\grad).';
                if norm(st) > pix, st = st*pix/norm(st); end
                p = p + st;
                if norm(p - p0) > 1.5*pix, p = p0 + (p - p0)*1.5*pix/norm(p - p0); break; end
                if norm(st) < 1e-9, break; end
            end
            est = [p; ic; wa; gau; par; xs(ix), zs(iz)];
            se(ie, :) = se(ie, :) + (est(:, 1).' - x0).^2;
            seZ(ie) = seZ(ie) + (p(2) - z0)^2;
            nOut(ie) = nOut(ie) + (abs(p(1) - x0) > fwhmX/2 || abs(p(2) - z0) > fwhmZ/2);
            if ip == 3 && t == nTrial && enrList(ie) == 20   % keep one example for the figure
                Yex = Y; xsEx = xs; zsEx = zs; x0Ex = x0; z0Ex = z0; estEx = est;
            end
        end
    end
    rmseX(ip, :, :) = reshape(sqrt(se/nTrial), [1, numel(enrList), nEst]);
    rmseZmle(ip, :) = sqrt(seZ.'/nTrial);
    outMLE(ip, :) = nOut.'/nTrial;
end
fprintf('PSF FWHM lateral %.2f lambda, axial %.2f lambda; f0 = %.1f MHz (lambda = %.0f um); %d trials\n', ...
    fwhmX, fwhmZ, f0/1e6, lam*1e6, nTrial);
fprintf('Lateral CRLB sigma_x = wx/sqrt(ENR); entries are RMSE / CRLB (lateral)\n');
for ip = 1:numel(pixList)
    fprintf('\npixel = lambda/%d (%.2f FWHM_x); peak-pixel SNR = ENR - %.1f dB\n', round(1/pixList(ip)), pixList(ip)/fwhmX, snrPkOff(ip));
    fprintf('  ENR(dB)  CRLB_x(lambda)  CRLB_x(um) ');
    fprintf('%7s', names{:}); fprintf('   MLE_z  out%%\n');
    for ie = 1:numel(enrList)
        fprintf('  %6d  %10.4f  %10.2f  ', enrList(ie), crlbX(ie), crlbX(ie)*lam*1e6);
        fprintf('%7.2f', squeeze(rmseX(ip, ie, :)).'/crlbX(ie));
        fprintf('  %6.2f  %4.1f\n', rmseZmle(ip, ie)/crlbZ(ie), 100*outMLE(ip, ie));
    end
end
fprintf('\nFWHM-lambda rule: ENR 20 dB -> sigma_x = lambda/%.1f; 30 dB -> lambda/%.1f\n', ...
    1/(wx/sqrt(100)), 1/(wx/sqrt(1000)));
fprintf('Elapsed %.1f s\n', toc);

%% Figures
figure('Position', [60 60 1150 360]);
subplot(1, 3, 1);
imagesc(xsEx, zsEx, 20*log10(abs(Yex)/max(abs(Yex(:))))); caxis([-30 0]); colormap(gray); axis image; colorbar;
hold on; plot(x0Ex, z0Ex, 'g+', 'MarkerSize', 12, 'LineWidth', 1.5);
plot(estEx(1,1), estEx(1,2), 'ro', estEx(4,1), estEx(4,2), 'cs', 'MarkerSize', 8);
xlabel('lateral x (\lambda)'); ylabel('axial z (\lambda)');
title('|IQ| dB, \lambda/8, ENR 20 dB: + true, o MLE, [] Gauss');
subplot(1, 3, 2);
ip = 2; sty = {'b-o', 'g-s', '-^', 'm-d', 'r-v', 'k:'};
for k = 1:nEst
    semilogy(enrList, squeeze(rmseX(ip, :, k)), sty{k}); hold on;
end
semilogy(enrList, crlbX, 'k-', 'LineWidth', 2);
xlabel('ENR (dB)'); ylabel('lateral RMSE (\lambda)'); grid on;
title('\lambda/4 pixels'); legend([names, {'CRLB'}], 'Location', 'northeast');
subplot(1, 3, 3);
ie = find(enrList == 30);
bar(squeeze(rmseX(:, ie, 1:5))/crlbX(ie));
set(gca, 'XTickLabel', {'\lambda/2', '\lambda/4', '\lambda/8'}); ylim([0 4]);
hold on; plot([0.5 3.5], [1 1], 'k--');
ylabel('RMSE / CRLB (lateral)'); title('ENR 30 dB'); legend(names(1:5), 'Location', 'northeast');

What to look for: the lateral CRLB column reads 13.08 µm at 20 dB and 4.14 µm at 30 dB (FWHM = λ at 5 MHz, λ/23.5 and λ/74.5). For ENR ≥ 20 dB the MLE sits at 0.95–1.03 times the bound at every pixel size, and the intensity centroid at 0.95–1.3. The 3-point Gaussian is 1.2–1.4 at λ/2 but about 2.2 at λ/4 and λ/8; the parabola reaches 5.7 at 40 dB with λ/2 pixels (bias floor), and the pixel-centre estimate grows as (pix/√12)/CRLB. At 16 dB, 1 % outliers already give 2–2.5 times the bound; at 10 dB about half the trials are outliers (out% column). The header line per pixel size gives the peak-pixel SNR offset (4.8 dB at λ/4 and λ/8).

m12_crlb_amplitude_mle.mDC level under Gaussian, Laplace and uniform noise; amplitude of a known pulse in real vs complex noise (the factor 2); bias of |â|² · change N, noise laws · runtime ≈ 2 s

%% m12_crlb_amplitude_mle.m
% Purpose : the Cramer-Rao bound in its simplest settings.
%   Part 1: DC level A in N noise samples. Sample mean, median and midrange
%           under Gaussian, Laplace and uniform noise of equal variance.
%           Prints N*var/sigma^2, so the Gaussian CRLB is 1 and the Laplace
%           CRLB is 0.5. Uniform noise violates the regularity condition
%           (its support moves with A): no CRLB, and the midrange beats 1/N.
%   Part 2: amplitude of a known pulse s in complex circular Gaussian noise
%           (IQ data). Real amplitude: var >= sigma^2/(2||s||^2).
%           Complex amplitude: E|a_hat - a|^2 >= sigma^2/||s||^2.
%           Same pulse in real noise of variance sigma^2: sigma^2/||s||^2.
%   Part 3: MLE invariance. The MLE of |a|^2 is |a_hat|^2, which is biased
%           upward by sigma^2/||s||^2.
% Change  : N, nTrial. Expect Laplace-median N*var to approach 0.5 slowly
%           (0.58 at N = 101).
% Expect  : N = 25: Gaussian mean 1.00, median ~1.54, midrange ~3.3;
%           Laplace mean 1.00, median ~0.67; uniform midrange ~0.214
%           (theory 2*3*N/((N+1)(N+2)) = 0.2137). Part 2 ratios ~1.00.
rng(1);

%% Part 1: DC level, three noise laws with unit variance
N = 25;            % samples per estimate
nTrial = 20000;    % Monte Carlo trials
A = 0;             % true level (any value; estimators are shift-equivariant)
lawName = {'Gaussian', 'Laplace ', 'uniform '};
crlbN = [1, 0.5, NaN];   % N*CRLB/sigma^2 (none for uniform)
res = zeros(3, 3);       % rows: law, cols: mean, median, midrange
estAll = cell(3, 1);
for law = 1:3
    switch law
        case 1
            w = randn(nTrial, N);
        case 2   % Laplace with variance 1: scale b = 1/sqrt(2)
            u = rand(nTrial, N) - 0.5;
            w = -(1/sqrt(2)) * sign(u) .* log(1 - 2*abs(u));
        case 3   % uniform on [-sqrt(3), sqrt(3)]: variance 1
            w = (rand(nTrial, N) - 0.5) * 2*sqrt(3);
    end
    y = A + w;
    est = [mean(y, 2), median(y, 2), (max(y, [], 2) + min(y, [], 2))/2];
    estAll{law} = est;
    res(law, :) = N * var(est - A);
end
fprintf('Part 1: N*var/sigma^2 for N = %d (%d trials)\n', N, nTrial);
fprintf('  noise      mean    median  midrange   N*CRLB\n');
for law = 1:3
    fprintf('  %s  %6.3f  %6.3f  %7.3f   %5.2f\n', lawName{law}, res(law, :), crlbN(law));
end
fprintf('  theory: median/Gaussian -> pi/2 = %.3f; Laplace median -> 0.5; uniform midrange = %.4f\n', ...
    pi/2, 6*N/((N+1)*(N+2)));

% N*var versus N for the Gaussian and Laplace cases
Ns = [3 5 9 17 33 65 129];
effN = zeros(numel(Ns), 4);  % Gauss mean, Gauss median, Laplace mean, Laplace median
for k = 1:numel(Ns)
    n = Ns(k); g = randn(5000, n);
    u = rand(5000, n) - 0.5; l = -(1/sqrt(2)) * sign(u) .* log(1 - 2*abs(u));
    effN(k, :) = n * [var(mean(g, 2)), var(median(g, 2)), var(mean(l, 2)), var(median(l, 2))];
end

%% Part 2: amplitude of a known complex pulse (IQ), real vs complex noise
fs = 40e6; f0 = 5e6; t = (-40:40)'/fs;          % time axis (s)
s = exp(-t.^2/(2*(0.15e-6)^2)) .* exp(1j*2*pi*f0*t);  % known complex pulse
sigma2 = 0.5;                                   % noise variance E|w|^2
a = 0.8;                                        % true real amplitude
M = 20000; E = sum(abs(s).^2);
W = sqrt(sigma2/2) * (randn(numel(s), M) + 1j*randn(numel(s), M));
Y = a*s + W;
proj = (s' * Y) / E;                            % s^H y / ||s||^2 for each trial
aReal = real(proj);                             % MLE when a is real (phase known)
aCplx = proj;                                   % MLE when a is complex
fprintf('\nPart 2: known pulse, ||s||^2 = %.2f, sigma^2 = %.2f\n', E, sigma2);
fprintf('  real a, complex noise : var = %.3e, CRLB sigma^2/(2||s||^2) = %.3e, ratio %.3f\n', ...
    var(aReal), sigma2/(2*E), var(aReal)/(sigma2/(2*E)));
fprintf('  complex a             : E|err|^2 = %.3e, CRLB sigma^2/||s||^2 = %.3e, ratio %.3f\n', ...
    mean(abs(aCplx - a).^2), sigma2/E, mean(abs(aCplx - a).^2)/(sigma2/E));
sr = real(s); Er = sum(sr.^2);                  % real pulse in real noise of variance sigma^2
yr = a*sr + sqrt(sigma2)*randn(numel(sr), M);
aR = (sr' * yr) / Er;
fprintf('  real pulse, real noise: var = %.3e, CRLB sigma^2/||s||^2 = %.3e, ratio %.3f\n', ...
    var(aR), sigma2/Er, var(aR)/(sigma2/Er));

%% Part 3: invariance is not unbiasedness
p = abs(aCplx).^2;                              % MLE of the power |a|^2
fprintf('\nPart 3: E[|a_hat|^2] = %.4f, |a|^2 = %.4f, bias = %.4f, predicted sigma^2/||s||^2 = %.4f\n', ...
    mean(p), a^2, mean(p) - a^2, sigma2/E);

%% Figures
figure('Position', [80 80 1000 380]);
subplot(1, 2, 1);
est = estAll{1}; edges = linspace(-1, 1, 61); ctr = (edges(1:end-1) + edges(2:end))/2;
h1 = histc(est(:, 1), edges); h2 = histc(est(:, 2), edges); h3 = histc(est(:, 3), edges);
bw = edges(2) - edges(1);
plot(ctr, h1(1:end-1)/(nTrial*bw), 'b-', ctr, h2(1:end-1)/(nTrial*bw), 'r-', ...
     ctr, h3(1:end-1)/(nTrial*bw), 'Color', [0.8 0.5 0], 'LineWidth', 1.4); hold on;
plot(ctr, sqrt(N/(2*pi)) * exp(-N*ctr.^2/2), 'k--', 'LineWidth', 1.2);
xlabel('estimate - A (noise std units)'); ylabel('density');
title(sprintf('Gaussian noise, N = %d', N));
legend('mean', 'median', 'midrange', 'CRLB pdf N(0,\sigma^2/N)');
subplot(1, 2, 2);
semilogx(Ns, effN(:, 1), 'b-o', Ns, effN(:, 2), 'r-o', Ns, effN(:, 3), 'b--s', Ns, effN(:, 4), 'r--s'); hold on;
semilogx(Ns, ones(size(Ns)), 'k:', Ns, 0.5*ones(size(Ns)), 'k-.');
xlabel('N (samples)'); ylabel('N \cdot var / \sigma^2'); ylim([0 2]);
title('Efficiency: Gaussian (solid) and Laplace (dashed) noise');
legend('mean, Gauss', 'median, Gauss', 'mean, Laplace', 'median, Laplace', 'CRLB Gauss', 'CRLB Laplace');

What to look for: \(N\cdot\operatorname{var}/\sigma^2\) at N = 25 of 1.00 (Gaussian mean), 1.55 (Gaussian median, \(\pi/2\) asymptotically), 0.68 (Laplace median, bound 0.5) and 0.215 for the uniform midrange (theory \(6N/((N+1)(N+2)) = 0.2137\)). Part 2 ratios to \(\sigma^2/(2\lVert s\rVert^2)\), \(\sigma^2/\lVert s\rVert^2\) and \(\sigma^2/\lVert s\rVert^2\) are 0.99, 1.00, 0.99. Part 3: \(\mathrm{E}\lvert\hat a\rvert^2\) exceeds \(\lvert a\rvert^2\) by 0.0474 against a predicted 0.0470.

m12_frequency_crlb.mComplex-tone MLE and lag-one estimator vs 6/(SNR·N(N²−1)), the threshold, and Kay's real-sinusoid bound · change N, the SNR grid · runtime ≈ 5 s

%% m12_frequency_crlb.m
% Purpose : frequency of a tone in noise (slow-time Doppler is this problem).
%   Part 1: complex tone y[n] = A exp(j(w n + phi)) + w[n], n = 0..N-1, in
%           circular complex Gaussian noise of variance sigma^2. The MLE is
%           the periodogram peak (zero-padded FFT, then a fine grid and a
%           parabola). Compared with the CRLB
%               var(w_hat) >= 6 / (SNR N (N^2 - 1)),  SNR = A^2/sigma^2,
%           and with the lag-one (Kasai) estimator. Shows the threshold.
%   Part 2: real sinusoid x[n] = A cos(2 pi f n + phi) + w[n] (Kay, Vol. I,
%           Example 3.14): var(f_hat) >= 12 / ((2 pi)^2 eta N (N^2 - 1)),
%           eta = A^2/(2 sigma^2). The exact MLE (projection on cos/sin)
%           reaches it; the plain periodogram peak does not at short N.
% Change  : N, the SNR grid, f0/PRF for the velocity conversion.
% Expect  : N = 32: MSE/CRLB ~1.0 for SNR >= -2 dB (N*SNR >= 13 dB); ~180x
%           at -4 dB; hundreds of times below. Kasai approaches
%           N(N+1)/(6(N-1)) = 5.7x the CRLB at high SNR and is far worse at
%           moderate SNR. Part 2: exact MLE ~1.0x; periodogram peak 1.4x,
%           4.3x, 35x at eta = 10, 20, 30 dB (bias from the -f image).
rng(2);

%% Part 1: complex tone, MLE and Kasai vs CRLB
N = 32;                       % samples (ensemble length)
snrdB = -10:2:24;             % per-sample SNR A^2/sigma^2 (dB)
nTrial = 1000;
L = 2^nextpow2(16*N);         % zero-padded FFT length
n = (0:N-1)';
nFine = 41;                   % fine grid points across +-1 FFT bin
mseMLE = zeros(size(snrdB)); mseKas = mseMLE; outFrac = mseMLE; crlb = mseMLE;
for k = 1:numel(snrdB)
    snr = 10^(snrdB(k)/10);
    w0 = (rand(1, nTrial) - 0.5) * 2*pi;          % true frequencies (rad/sample)
    ph = rand(1, nTrial) * 2*pi;
    Y = exp(1j*(n*w0 + ones(N, 1)*ph)) + sqrt(1/(2*snr)) * (randn(N, nTrial) + 1j*randn(N, nTrial));
    P = abs(fft(Y, L)).^2;
    [~, kmax] = max(P);
    wc = 2*pi*(kmax - 1)/L;                       % coarse peak
    errM = zeros(1, nTrial);
    for t = 1:nTrial
        wg = wc(t) + linspace(-2*pi/L, 2*pi/L, nFine);   % fine grid
        J = abs(exp(-1j*n*wg).' * Y(:, t)).^2;            % periodogram on the grid
        [~, j] = max(J); j = min(max(j, 2), nFine - 1);
        dlt = 0.5*(J(j-1) - J(j+1)) / (J(j-1) - 2*J(j) + J(j+1));
        wHat = wg(j) + dlt*(wg(2) - wg(1));
        errM(t) = angle(exp(1j*(wHat - w0(t))));  % wrap to (-pi, pi]
    end
    errK = angle(sum(Y(2:end, :) .* conj(Y(1:end-1, :)), 1) .* exp(-1j*w0));
    crlb(k) = 6 / (snr * N * (N^2 - 1));
    mseMLE(k) = mean(errM.^2); mseKas(k) = mean(errK.^2);
    outFrac(k) = mean(abs(errM) > 0.5*2*pi/N);    % beyond half a mainlobe
end
fprintf('Part 1: complex tone, N = %d, %d trials per SNR\n', N, nTrial);
fprintf('  SNR(dB)  N*SNR(dB)  MLE/CRLB  Kasai/CRLB  outliers\n');
for k = 1:numel(snrdB)
    fprintf('  %6.0f  %8.1f  %9.2f  %9.2f  %8.3f\n', snrdB(k), snrdB(k) + 10*log10(N), ...
        mseMLE(k)/crlb(k), mseKas(k)/crlb(k), outFrac(k));
end
% Doppler conversion at one operating point
c = 1540; f0 = 5e6; PRF = 5e3; M = 16; snr = 10;   % 10 dB per pulse
sw = sqrt(6/(snr*M*(M^2 - 1)));
fprintf('  Doppler example: M = %d, SNR = 10 dB: sigma_w = %.4f rad/sample, sigma_fD = %.2f Hz, sigma_v = %.2f mm/s\n', ...
    M, sw, sw*PRF/(2*pi), c*sw*PRF/(2*pi)/(2*f0)*1e3);

%% Part 2: real sinusoid, exact MLE vs periodogram peak vs Kay's CRLB
N2 = 32; n2 = (0:N2-1)'; etadB = [10 20 30]; nT2 = 2000; L2 = 512;
fprintf('\nPart 2: real sinusoid, N = %d, f0 ~ U(0.1, 0.4) cycles/sample\n', N2);
for k = 1:numel(etadB)
    eta = 10^(etadB(k)/10); s2 = 1/(2*eta);              % A = 1
    f0t = 0.1 + 0.3*rand(1, nT2); ph = 2*pi*rand(1, nT2);
    X = cos(2*pi*n2*f0t + ones(N2, 1)*ph) + sqrt(s2)*randn(N2, nT2);
    eP = zeros(1, nT2); eX = eP;
    for t = 1:nT2
        x = X(:, t);
        Pq = abs(fft(x, L2)).^2; [~, kk] = max(Pq(2:L2/2)); kk = kk + 1;
        fg = (kk - 1)/L2 + linspace(-1/L2, 1/L2, nFine);
        E = exp(-1j*2*pi*n2*fg);
        Jp = abs(E.' * x).^2;                           % periodogram on the grid
        C = real(E); S = -imag(E);                       % cos and sin columns
        xc = C.'*x; xs = S.'*x; cc = sum(C.^2).'; ss = sum(S.^2).'; cs = sum(C.*S).';
        Jx = (ss.*xc.^2 - 2*cs.*xc.*xs + cc.*xs.^2) ./ (cc.*ss - cs.^2);   % ||P_H x||^2
        [~, j] = max(Jp); j = min(max(j, 2), nFine - 1);
        eP(t) = fg(j) + 0.5*(Jp(j-1) - Jp(j+1))/(Jp(j-1) - 2*Jp(j) + Jp(j+1))*(fg(2) - fg(1)) - f0t(t);
        [~, j] = max(Jx); j = min(max(j, 2), nFine - 1);
        eX(t) = fg(j) + 0.5*(Jx(j-1) - Jx(j+1))/(Jx(j-1) - 2*Jx(j) + Jx(j+1))*(fg(2) - fg(1)) - f0t(t);
    end
    crlbR = 12 / ((2*pi)^2 * eta * N2 * (N2^2 - 1));
    fprintf('  eta = %2d dB: exact MLE var/CRLB = %.2f, periodogram peak MSE/CRLB = %.2f\n', ...
        etadB(k), var(eX)/crlbR, mean(eP.^2)/crlbR);
end

%% Figure
figure('Position', [80 80 900 380]);
subplot(1, 2, 1);
semilogy(snrdB, sqrt(mseMLE), 'bo-', snrdB, sqrt(mseKas), 's-', 'Color', [0.8 0.5 0]); hold on;
semilogy(snrdB, sqrt(crlb), 'k-', 'LineWidth', 1.4);
semilogy(snrdB, pi/sqrt(3)*ones(size(snrdB)), 'k:');
xlabel('per-sample SNR (dB)'); ylabel('RMS error (rad/sample)');
title(sprintf('Complex tone, N = %d', N)); legend('periodogram MLE', 'Kasai (lag one)', 'CRLB', 'uniform guess pi/sqrt(3)');
subplot(1, 2, 2);
plot(snrdB + 10*log10(N), outFrac, 'bo-'); grid on;
xlabel('integrated SNR N\cdotSNR (dB)'); ylabel('outlier fraction');
title('Threshold: outliers beyond half a mainlobe');

What to look for: for N = 32 the MLE/CRLB ratio is 0.94–1.09 from −2 dB per sample (13.1 dB integrated) up, 182 at −4 dB and above 500 below; the outlier fraction jumps from 0 to 2 % to 12 % across 13 → 11 → 9 dB integrated. Lag-one falls toward 5.7 (5.7–6.3 at 22–24 dB). Doppler example: 9.65 Hz and 1.49 mm/s. Part 2: the exact real-sinusoid MLE sits at 0.99–1.04 of Kay's bound while the plain periodogram peak is 1.4, 4.3 and 35 times above it at η = 10, 20, 30 dB (image interference at N = 32).

m12_nuisance_fisher.mStd inflation from unknown amplitude, width, background, a neighbour bubble, carrier phase; CRLB on sound speed with z, t₀, x_s unknown · change offsets, separations, array, σ_τ · runtime ≈ 2 s

%% m12_nuisance_fisher.m
% Purpose : what unknown (nuisance) parameters cost, from the Fisher matrix.
%   The bound on parameter 1 with nuisances is [inv(I)](1,1) >= 1/I(1,1).
%   We print the std inflation  sqrt([inv(I)](1,1) * I(1,1))  (1 = free).
%   Noise: circular complex Gaussian, E|w|^2 = 1 per sample, so
%   I = 2 Re(J' * J) with J the Jacobian of the noise-free signal.
%   Part 1: one bubble, Gaussian PSF (FWHM = 1): unknown complex amplitude,
%           width, complex background; centred vs near the window edge.
%   Part 2: a second bubble at separation d (FWHM units), relative phase psi,
%           its position and amplitude unknown. Inflation vs d.
%   Part 3: axial localisation of a carrier-modulated echo: known vs unknown
%           phase. With unknown phase the carrier term cancels exactly.
%   Part 4: sound speed from receive delays of one point target
%           (0 deg plane-wave transmit, 128 elements, pitch 0.3 mm,
%           1 ns independent delay jitter per element).
% Change  : window half-width, offsets, phases, the array, sigma_tau.
% Expect  : Part 1: 1.000 for a centred bubble with every nuisance unknown;
%           ~4.2 when the bubble sits at the edge of a +-1.5 FWHM window.
%           Part 2: in phase 3.96 / 1.39 / 1.00 at d = 0.5 / 1 / 2 FWHM,
%           quadrature 1.81 / 1.11 / 1.00; same values for a weaker neighbour.
%           Part 3: unknown-phase Fisher = envelope-only Fisher (ratio 1.000).
%           Part 4 at z = 30 mm: sigma_c = 0.0034 (z known), 0.065 (z
%           unknown), 2.83 m/s (z and t0 unknown).
rng(1);
fw = 1; w = fw/(2*sqrt(2*log(2))); pix = 1/16;
gfun = @(x, x0) exp(-(x - x0).^2/(2*w^2));
first = @(v) v(1);
inv11 = @(M) first(M \ [1; zeros(size(M, 1) - 1, 1)]);     % [inv(M)](1,1)
infl = @(J) sqrt(inv11(2*real(J'*J)) * 2*real(J(:,1)'*J(:,1)));

%% Part 1: single bubble
fprintf('Part 1: one bubble, std inflation of x1 (1.000 = nuisance is free)\n');
cases = {'centred, window +-3 FWHM', 3, 0; 'offset 1.0, window +-1.5', 1.5, 1.0; 'offset 1.5, window +-1.5 (edge)', 1.5, 1.5};
for k = 1:size(cases, 1)
    win = cases{k, 2}; off = cases{k, 3};
    x = (-round(win/pix):round(win/pix))'*pix; g = gfun(x, off);
    dx = g.*(x - off)/w^2; dw = g.*(x - off).^2/w^3; one = ones(size(x));
    J1 = [dx, g, 1j*g];                       % x1, Re a, Im a
    J2 = [J1, dw];                            % + width
    J3 = [J2, one, 1j*one];                   % + complex background
    fprintf('  %-34s amp: %.3f   amp+width: %.3f   amp+width+bg: %.3f\n', cases{k, 1}, infl(J1), infl(J2), infl(J3));
end

%% Part 2: a neighbour bubble
dList = [0.25 0.5 0.75 1 1.5 2 3];
psiList = [0 90 180]; relList = [1 0.3];
x = (-4/pix:4/pix)'*pix; ratio2 = zeros(numel(psiList), numel(dList));
fprintf('\nPart 2: second bubble at separation d (x2, a2 unknown; a1 unknown)\n');
fprintf('  psi(deg) |a2|  '); fprintf('  d=%-5.2f', dList); fprintf('\n');
for ir = 1:numel(relList)
    for ip = 1:numel(psiList)
        a2 = relList(ir)*exp(1j*psiList(ip)*pi/180);
        r = zeros(size(dList));
        for k = 1:numel(dList)
            g1 = gfun(x, 0); g2 = gfun(x, dList(k));
            J = [g1.*x/w^2, g1, 1j*g1, a2*g2.*(x - dList(k))/w^2, g2, 1j*g2];
            r(k) = infl(J);
        end
        if ir == 1, ratio2(ip, :) = r; end
        fprintf('  %6d   %4.1f ', psiList(ip), relList(ir)); fprintf('  %7.2f', r); fprintf('\n');
    end
end
dS = [0.05 0.1 0.2]; rS = zeros(2, 3);
for k = 1:3
    for ip = 1:2
        a2 = exp(1j*(ip - 1)*pi/2); g1 = gfun(x, 0); g2 = gfun(x, dS(k));
        rS(ip, k) = infl([g1.*x/w^2, g1, 1j*g1, a2*g2.*(x - dS(k))/w^2, g2, 1j*g2]);
    end
end
fprintf('  small-d slope of log(inflation) vs log(d): in phase %.2f, quadrature %.2f\n', ...
    log(rS(1,3)/rS(1,1))/log(4), log(rS(2,3)/rS(2,1))/log(4));

%% Part 3: carrier with known vs unknown phase (axial, z in wavelengths)
fwz = 0.75; wz = fwz/(2*sqrt(2*log(2))); z = (-3:1/32:3)'; k0 = 2*pi*2;   % two-way carrier: 2/lambda
g = exp(-z.^2/(2*wz^2)); s = g.*exp(1j*k0*z);
ds = (-z/wz^2 + 1j*k0).*s;                    % d s / d z0 (sign irrelevant for Fisher)
Iknown = 2*real(ds'*ds);                      % amplitude and phase known
Junk = [ds, s, 1j*s];                         % complex amplitude unknown
Ieff = 1/inv11(2*real(Junk'*Junk));
dg = -z/wz^2.*g; Ienv = 2*real(dg'*dg);       % envelope only, no carrier
fprintf('\nPart 3: axial Fisher (FWHM %.2f lambda, carrier 2/lambda)\n', fwz);
fprintf('  phase known: %.1f   phase unknown: %.2f   envelope only: %.2f   ratio unknown/envelope = %.4f\n', ...
    Iknown, Ieff, Ienv, Ieff/Ienv);
fprintf('  std gain from a known phase: %.1fx\n', sqrt(Iknown/Ieff));

%% Part 4: sound speed from receive delays of one point target
c = 1540; N = 128; pitch = 0.3e-3; sig = 1e-9; D = N*pitch;
xe = ((0:N-1)' - (N-1)/2)*pitch;
zList = [10 20 30 45 60 80]*1e-3;
sets = {'c only', 'c, z', 'c, t0', 'c, z, t0', 'c, z, xs, t0'};
sigc = zeros(numel(zList), numel(sets));
for iz = 1:numel(zList)
    zz = zList(iz); r = sqrt(zz^2 + xe.^2);
    dc = -(zz + r)/c^2; dz = (1 + zz./r)/c; dxs = xe./(r*c); dt0 = ones(N, 1);
    Js = {dc, [dc dz], [dc dt0], [dc dz dt0], [dc dz dxs dt0]};
    for k = 1:numel(Js)
        Ci = inv(Js{k}'*Js{k}/sig^2); sigc(iz, k) = sqrt(Ci(1, 1));
    end
end
fprintf('\nPart 4: sigma_c (m/s), %d elements, pitch %.1f mm, sigma_tau = %.1f ns\n', N, pitch*1e3, sig*1e9);
fprintf('  z(mm)  F#  '); fprintf('%14s', sets{:}); fprintf('   approx(c,z)\n');
for iz = 1:numel(zList)
    fprintf('  %4.0f  %4.2f', zList(iz)*1e3, zList(iz)/D); fprintf('%14.4f', sigc(iz, :));
    fprintf('   %9.4f\n', sqrt(180)*sig*zList(iz)*c^2/(sqrt(N)*D^2));
end

%% Figures
figure('Position', [80 80 950 360]);
subplot(1, 2, 1);
semilogy(dList, ratio2(1, :), 'b-o', dList, ratio2(2, :), 'r-s'); hold on; semilogy(dList, ones(size(dList)), 'k:');
xlabel('separation d (FWHM)'); ylabel('std inflation of x_1'); grid on;
title('Neighbour bubble, position and amplitude unknown'); legend('in phase', 'quadrature');
subplot(1, 2, 2);
semilogy(zList*1e3, sigc(:, 1), 'k-o', zList*1e3, sigc(:, 2), 'b-s', zList*1e3, sigc(:, 3), 'm-^', zList*1e3, sigc(:, 4), 'r-d');
xlabel('depth z (mm)'); ylabel('\sigma_c (m/s)'); grid on;
title('Sound speed from one target, 1 ns jitter'); legend(sets(1:4), 'Location', 'northwest');

What to look for: Part 1: 1.000 for a centred bubble with every nuisance unknown, 1.55 / 4.12 / 4.20 at the window edge. Part 2: 14.55, 3.96, 2.03, 1.38, 1.04 at d = 0.25…1.5 FWHM in phase, 3.45, 1.81, 1.31, 1.11, 1.01 in quadrature, identical for \(\lvert a_2\rvert\) = 0.3; small-d slopes −1.99 and −0.99. Part 3: unknown-phase Fisher equals envelope-only Fisher (ratio 1.0000), known phase 5.7× better in σ. Part 4 at 30 mm: \(\sigma_c\) = 0.0034, 0.065, 0.124 and 2.83 m/s.

m12_shrinkage_mse.mScalar shrinkage MSE vs the CRLB, exact and Monte Carlo James–Stein risk, posterior mean = Bayesian CRB · change c, p, τ · runtime ≈ 3 s

%% m12_shrinkage_mse.m
% Purpose : biased estimators can have lower MSE than the CRLB, and the
%           Bayesian view explains when that is legitimate.
%   Part 1: scalar shrinkage theta_hat = c * ybar, ybar ~ N(theta, s^2),
%           s^2 = sigma^2/N (the CRLB). MSE = c^2 s^2 + (1-c)^2 theta^2,
%           below the CRLB only for |theta| < s*sqrt((1+c)/(1-c)).
%   Part 2: James-Stein for p >= 3 means, y ~ N(theta, I_p):
%           theta_JS = (1 - (p-2)/||y||^2) y. Exact risk
%           R = p - (p-2)^2 E[1/(p-2+2K)], with K the Poisson mixing index
%           (mean ||theta||^2/2) of the noncentral chi-square ||y||^2,
%           always < p (the MLE risk). Positive-part version is better still.
%   Part 3: Bayes. Prior theta ~ N(0, tau^2): the posterior mean is
%           shrinkage with c = tau^2/(tau^2 + s^2); its Bayesian MSE equals
%           the Bayesian CRB 1/(1/s^2 + 1/tau^2) (data + prior information).
% Change  : c, p, tau, the theta grid.
% Expect  : Part 1, c = 0.8: break-even |theta| = 3.00 s. Part 2: risk 2.00
%           at theta = 0 for every p; MC matches the exact formula to ~1 %.
%           Part 3: Bayesian MSE = Bayesian CRB (ratio 1.00).
rng(1);

%% Part 1: scalar shrinkage
s = 1;                                   % CRLB std of the unbiased estimate
cList = [0.5 0.8 0.95];
th = linspace(-5, 5, 401);
fprintf('Part 1: scalar shrinkage c*ybar, MSE in units of the CRLB\n');
for c = cList
    fprintf('  c = %.2f: MSE(0) = %.3f, break-even |theta| = %.2f s, MSE(3 s) = %.2f\n', ...
        c, c^2, sqrt((1 + c)/(1 - c)), c^2 + (1 - c)^2*9);
end
M = 200000; theta0 = 1.0; c = 0.8;
yb = theta0 + s*randn(M, 1);
fprintf('  Monte Carlo, theta = %.1f s, c = %.1f: MSE = %.4f (theory %.4f)\n', theta0, c, ...
    mean((c*yb - theta0).^2), c^2 + (1 - c)^2*theta0^2);

%% Part 2: James-Stein
jsRisk = @(p, lam) p - (p - 2)^2 * sum(exp(-lam/2 + (0:400)*log(max(lam/2, realmin)) - gammaln(1:401)) ./ (p - 2 + 2*(0:400)));
pList = [3 10 50]; normList = [0 1 2 3];   % ||theta|| / sqrt(p)
fprintf('\nPart 2: James-Stein risk (MLE risk = p)\n');
fprintf('     p  ||th||/sqrt(p)   exact     MC JS   MC JS+\n');
riskCurve = zeros(numel(pList), 31); nn = linspace(0, 3, 31);
for ip = 1:numel(pList)
    p = pList(ip);
    for k = 1:numel(normList)
        theta = zeros(1, p); theta(1) = normList(k)*sqrt(p);
        Y = theta + randn(40000, p);           % implicit expansion (R2016b+, Octave)
        n2 = sum(Y.^2, 2);
        JS = (1 - (p - 2)./n2) .* Y; JSp = max(0, 1 - (p - 2)./n2) .* Y;
        fprintf('  %4d  %10.1f  %9.3f %9.3f %8.3f\n', p, normList(k), jsRisk(p, sum(theta.^2)), ...
            mean(sum((JS - theta).^2, 2)), mean(sum((JSp - theta).^2, 2)));
    end
    for k = 1:numel(nn), riskCurve(ip, k) = jsRisk(p, nn(k)^2*p)/p; end
end

%% Part 3: Bayes posterior mean = shrinkage; Bayesian CRB
tau = 1.5; cB = tau^2/(tau^2 + s^2);
thetaB = tau*randn(M, 1); yB = thetaB + s*randn(M, 1);
mseB = mean((cB*yB - thetaB).^2); bcrb = 1/(1/s^2 + 1/tau^2);
fprintf('\nPart 3: prior std tau = %.1f s: c = %.3f, Bayesian MSE = %.4f, Bayesian CRB = %.4f, ratio %.3f\n', ...
    tau, cB, mseB, bcrb, mseB/bcrb);
fprintf('  same estimator, MSE at theta = 3 tau = %.2f (CRLB 1): the prior must be right\n', ...
    cB^2 + (1 - cB)^2*(3*tau)^2);

%% Figures
figure('Position', [80 80 950 360]);
subplot(1, 2, 1);
for c = cList, plot(th, c^2 + (1 - c)^2*th.^2, 'LineWidth', 1.3); hold on; end
plot(th, ones(size(th)), '--', 'Color', [0.6 0 0], 'LineWidth', 2); ylim([0 3]); grid on;
xlabel('\theta / s  (s = CRLB std)'); ylabel('MSE / CRLB');
title('Shrinkage c\cdot y: below the CRLB near \theta = 0 only');
legend('c = 0.5', 'c = 0.8', 'c = 0.95', 'unbiased (CRLB)', 'Location', 'north');
subplot(1, 2, 2);
plot(nn, riskCurve.', 'LineWidth', 1.3); hold on; plot(nn, ones(size(nn)), '--', 'Color', [0.6 0 0], 'LineWidth', 2);
xlabel('||\theta|| / sqrt(p)'); ylabel('risk / p'); ylim([0 1.1]); grid on;
title('James-Stein dominates the MLE for p >= 3');
legend('p = 3', 'p = 10', 'p = 50', 'MLE', 'Location', 'southeast');

What to look for: c = 0.8 breaks even at |θ| = 3.00 s; Monte Carlo MSE 0.6796 vs 0.6800. James–Stein risk is 2.000 at θ = 0 for p = 3, 10 and 50 and always below p; the positive-part version is lower still (1.61, 1.24, 1.12 at θ = 0). Part 3: Bayesian MSE 0.6925 against the Bayesian CRB 0.6923, and 2.40 times the CRLB for the same estimator at θ = 3τ.

12.13Self-check

Q1A tracker has 2 µm bias and 3 µm standard deviation. What is its RMSE, and what happens to it with 4× more independent frames averaged?

\(\mathrm{RMSE} = \sqrt{2^2 + 3^2} = 3.61\) µm. Averaging 4 frames halves the std to 1.5 µm but leaves the bias: \(\sqrt{4 + 2.25} = 2.5\) µm. The bias becomes the floor.

Q2A bubble PSF has a lateral FWHM of 200 µm and ENR = 25 dB. What is the lateral CRLB?

\(w = 200/2.355 = 84.9\) µm; \(\sqrt{\mathrm{ENR}} = \sqrt{316.2} = 17.8\); \(\sigma_x \ge 4.8\) µm (FWHM/41.9).

Q3You resample the same beamformed IQ image from λ/4 to λ/16 pixels. Does the localization CRLB change? Does your localizer's error?

The bound does not change: noise density and PSF are the same, and the extra samples carry correlated noise. A matched-filter MLE or intensity centroid stays near the bound; a 3-point fit gets worse (it now sees a smaller part of the PSF). Only a simulation that adds independent noise per pixel would report an improvement.

Q4Slow-time ensemble of M = 16 at 10 dB per sample, PRF 5 kHz, \(f_0\) = 5 MHz. Lower bound on the velocity std?

\(\sigma_\omega = \sqrt{6/(10\cdot 16\cdot 255)} = 0.0121\) rad/sample → \(\sigma_f = 0.0121\cdot 5000/(2\pi) = 9.65\) Hz → \(\sigma_v = 1540\cdot 9.65/(2\cdot 5\times10^6) = 1.49\) mm/s. Doubling M to 32 divides it by 2.83.

Q5Why does an unknown complex amplitude cost nothing for the lateral position of a centred, symmetric PSF?

The Fisher cross term is \(\propto \mathrm{Re}\sum(a\,g')^*\,g\): \(g'\) is odd, \(g\) is even, so the sum is zero. The position derivative is orthogonal to the amplitude derivatives, and \([\mathbf{I}^{-1}]_{11} = 1/I_{11}\). Off-centre in a finite window, or with an asymmetric PSF, the symmetry breaks and the cost appears.

Q6A frequency estimator's RMSE is on the CRLB at 0 dB per sample but 100× above it at −4 dB, with N = 32. What is happening, and which number predicts it?

The threshold effect: noise peaks in other bins beat the signal bin and the estimate jumps (outliers), which the local CRLB cannot describe. The relevant number is the integrated SNR \(N\cdot\mathrm{SNR}\) (15 dB vs 11 dB here); global bounds such as Ziv–Zakai capture it.

Q7Shrinkage with c = 0.5: below the CRLB for which true values, and which prior makes it optimal?

\(\mathrm{MSE} = 0.25s^2 + 0.25\theta^2 < s^2\) for \(\lvert\theta\rvert < \sqrt{3}\,s = 1.73s\). It is the posterior mean for a prior std \(\tau = s\sqrt{c/(1 - c)} = s\), with Bayes MSE \(c\,s^2 = 0.5s^2\).

12.14Figures from the literature

Published figures with real data that show the quantities of this module: localization precision against the PSF, the pixel versus the precision, the cost of neighbours, and how resolution of the final map is measured. Links only; figure numbers checked against the publisher's pages.

Ultrafast ULM in the rat brain: precision and resolution measured in vivo

C. Errico, J. Pierre, S. Pezet, et al., M. Tanter, "Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging," Nature 527, 499–502, 2015. Fig. 2 ("Spatial resolution and quantification of uULM in the rat brain cortex through a thinned skull window"). doi:10.1038/nature16066

Look at how the vessel profiles on the ULM map compare with the diffraction-limited PSF of the probe, and at the gap between the precision of single localizations and the width of the reconstructed vessels. The first number is the \(w/\sqrt{\mathrm{ENR}}\) of (12.11); the second adds tracking, motion and finite counts.

A ULM map reconstructed on 5 µm pixels

V. Hingot, C. Errico, B. Heiles, L. Rahal, M. Tanter, O. Couture, "Microvascular flow dictates the compromise between spatial resolution and acquisition time in Ultrasound Localization Microscopy," Scientific Reports 9, 2456, 2019 (open access, CC BY 4.0). Fig. 2 (coronal slice of the rat brain, 5 µm pixels, 240 s). doi:10.1038/s41598-018-38349-x

The display pixel is far below λ, but each pixel is a histogram of localizations: its value is meaningful only once enough independent localizations fall in it. Read it with Section 12.5 in mind: the pixel neither sets nor improves the bound; it sets how many detections a vessel needs before its texture stops being jitter.

How long you must acquire before small vessels saturate

Hingot et al., Sci. Rep. 9, 2456, 2019 (as above). Fig. 5 ("Determination of the influence of the acquisition time on vessel reconstruction"). doi:10.1038/s41598-018-38349-x

How the reconstruction of individual vessels fills in as acquisition time grows. The CRLB per localization is only the first factor in map quality; the number of localizations per resolution cell, set by bubble concentration and flow, is the second, and it is the one that costs minutes.

Measuring the resolution of a human transcranial ULM map

C. Demené, J. Robin, A. Dizeux, et al., F. Perren, "Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients," Nature Biomedical Engineering 5, 219–228, 2021. Fig. 2 ("Estimation of the ULM resolution in a standard imaging case"). doi:10.1038/s41551-021-00697-x

An example of reporting the resolution of the final map rather than localization precision, here through the skull. Compare the quoted resolution with the probe's FWHM and with \(w/\sqrt{\mathrm{ENR}}\): the ratio tells you how much of the map's quality is limited by precision and how much by counts, motion and aberration.

Localization error and missed bubbles versus concentration

Y. Shin, M. R. Lowerison, Y. Wang, et al., P. Song, "Context-aware deep learning enables high-efficacy localization of high concentration microbubbles for super-resolution ultrasound localization microscopy," Nature Communications 15, 2932, 2024 (open access). Fig. 2b (detection accuracy, missing rate and localization error against microbubble concentration, simulation). doi:10.1038/s41467-024-47154-2

As concentration rises, PSFs overlap and conventional localization loses bubbles and precision. This is the nuisance-parameter cost of Section 12.6 (inflation 3.96 at 0.5 FWHM in phase) plus the detection failure the Fisher matrix does not show. A learned localizer moves the curve, which is the Bayesian-prior effect of Section 12.7.

12.15Go deeper

S. M. Kay, Fundamentals of Statistical Signal Processing, Vol. I: Estimation Theory, Prentice Hall, 1993. Chapters 3 (CRLB, incl. Example 3.14), 7 (MLE), 8 (least squares), 10–11 (Bayesian).
H. L. Van Trees, Detection, Estimation, and Modulation Theory, Part I, Wiley, 1968 (2nd ed. with K. L. Bell and Z. Tian, 2013). Classical and Bayesian bounds.
H. L. Van Trees and K. L. Bell (eds.), Bayesian Bounds for Parameter Estimation and Nonlinear Filtering/Tracking, Wiley–IEEE Press, 2007. Ziv–Zakai, Weiss–Weinstein and Bayesian CRB papers in one place.
D. C. Rife and R. R. Boorstyn, "Single-tone parameter estimation from discrete-time observations," IEEE Trans. Information Theory, 1974. Frequency CRLB, periodogram MLE and its threshold.
J. Ziv and M. Zakai, "Some lower bounds on signal parameter estimation," IEEE Trans. Information Theory, 1969.
C. Stein, "Inadmissibility of the usual estimator for the mean of a multivariate normal distribution," Proc. Third Berkeley Symposium, 1956; W. James and C. Stein, "Estimation with quadratic loss," Proc. Fourth Berkeley Symposium, 1961.
C. Errico et al., "Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging," Nature, 2015.
B. Heiles et al., "Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy," Nature Biomedical Engineering, 2022 (the PALA benchmark).
K. Christensen-Jeffries et al., "Super-resolution ultrasound imaging," Ultrasound in Medicine and Biology, 2020. Review of ULM methods and their error sources.
W. F. Walker and G. E. Trahey, "A fundamental limit on delay estimation using partially correlated speckle signals," IEEE TUFFC, 1995 (the delay bound used in Module 11 and as σ_τ in Demo 12.5).
C. Kasai, K. Namekawa, A. Koyano and R. Omoto, "Real-time two-dimensional blood flow imaging using an autocorrelation technique," IEEE Trans. Sonics and Ultrasonics, 1985.

← 11 Correlation & delay 13 Inverse problems →