Report an errorErrata
Part IV · Module 27

Spectral estimation & quantitative ultrasound

Every QUS number (attenuation coefficient, backscatter coefficient, effective scatterer diameter, spectral slope) is a function of an estimated power spectrum of a short, speckled RF gate. This module rebuilds the estimators, states their bias and variance in closed form, and carries the units (Np vs dB, amplitude vs power, one-way vs two-way) through every step to the attenuation and size estimates.

You will be able to

  • Predict the bias (window convolution) and variance (\(\chi^2_2\), 5.57 dB per bin) of a periodogram, and what Welch averaging, multitapers and AR models change
  • Attach a degrees-of-freedom count and a \(\chi^2\) confidence interval to any averaged spectrum
  • Derive the \(-4\alpha(f)\Delta z\) log spectral ratio (power, two-way, \(\alpha\) in Np) and convert it to dB/cm/MHz without a factor-2 or 8.686 error
  • Choose between the spectral log difference, the reference phantom regression and the centroid downshift, and predict each one's bias sources and standard deviation from \(K\), \(B\), \(T\) and the ROI length
  • Fit a Gaussian form factor for the effective scatterer diameter, know the \(ka\) range where it is identifiable, and see how an attenuation error leaks into it

Assumes

  • Module 04 · Sampling, Fourier & LTI systems: DFT, windows, zero padding, leakage
  • Module 07 · Random vectors, covariance & correlation: Wiener–Khinchin (the speckle power spectrum is the system spectrum)
  • Module 08 · Analytic signal & IQ demodulation: attenuation downshift of a Gaussian spectrum, Np vs dB
  • Module 12 · Estimation theory & the Cramér–Rao bound: bias, variance, CRLB vocabulary
Exact formulaclosed form plotted as isSimulationnumerical solution run in this page (JS) or in the Octave scriptsApprox. modelsimplified physics; assumptions under the badgeSchematicdrawing for intuition, not computed, not to scaleSpectra here come from three codes: FFT periodograms of simulated Gaussian RF (white noise through the stated transfer function), the multiplicative speckle model (an expected spectrum times the periodogram of a white-noise gate, averaged over K lines), and a toy focal-gain model for diffraction. No Field II or k-Wave output is shown.

27.1The periodogram: biased by the window, never consistent

Take \(N\) samples \(x[n]\) of a zero-mean stationary process at sampling rate \(f_s\). Its power spectral density is the DTFT of the autocorrelation (Module 07). Use per-sample units so that unit-variance white noise has \(S(f) = 1\); divide by \(f_s\) for per-hertz units.

\[ S(f) = \sum_{l} r[l]\, e^{-j 2\pi f l / f_s}, \qquad r[l] = \mathrm{E}\{ x[n]\, x[n+l] \} \]
(27.1)
\(x[n]\)
zero-mean stationary samples, \(n = 0 \dots N-1\)
\(f_s\)
sampling rate, MHz (40 MHz in the demos)
\(r[l]\)
autocorrelation at lag \(l\) samples
\(S(f)\)
power spectral density, per-sample units
\(l\)
lag, samples
\(\mathrm{E}\{\cdot\}\)
expected value
\(j\)
imaginary unit
\(\sum_l\)
sum over all lags
\(n\)
sample index
\[ P_w(f) = \frac{\left\lvert \sum_{n=0}^{N-1} w[n]\, x[n]\, e^{-j2\pi f n/f_s} \right\rvert^2}{\sum_n w[n]^2} \]
(27.2)
\(w[n]\)
data window (taper), dimensionless
\(P_w(f)\)
windowed periodogram, per-sample units
\(x[n],\ f_s,\ N\)
as in (27.1); N is the number of samples in the window
\(\lvert\cdot\rvert^2\)
squared magnitude
\(\sum_n w[n]^2\)
window energy (normalization)
The windowed periodogram. Rectangular \(w = 1\) gives the classical \((1/N)\lvert X(f)\rvert^2\). Normalizing by \(\sum w^2\) makes \(\mathrm{E}\{P_w\} = S\) for white noise for any window.

Bias. Multiplying by \(w\) in time convolves in frequency, so the expected periodogram is the true spectrum smeared by the window's power response (unit area):

\[ \mathrm{E}\{ P_w(f) \} = \int S(\nu)\, H_w(f - \nu)\, d\nu, \qquad H_w(f) = \frac{\lvert W(f)\rvert^2}{f_s \sum w^2} \]
(27.3)
\(W(f)\)
DTFT of the window
\(H_w(f)\)
spectral window (power response, unit area), 1/MHz
\(S(\nu)\)
true spectrum, integration variable \(\nu\) (MHz), as in (27.1)
\(P_w,\ w\)
as in (27.2)
\(\mathrm{E}\{\cdot\}\)
expected value
\(f,\ \nu\)
analysis frequency and integration frequency, MHz
\(f_s\)
sampling rate, MHz
Rectangular window: \(H\) is the Fejér kernel, mainlobe \(\pm f_s/N\), first sidelobe −13.3 dB, sidelobes decaying 6 dB/octave. Two consequences: resolution about \(f_s/N\), and leakage that fills in the low-power parts of a steep spectrum. Hann: sidelobes −31.5 dB, decaying 18 dB/octave, mainlobe twice as wide (ENBW \(1.5 f_s/N\)).

Variance. For a Gaussian process, away from 0 and \(f_s/2\) and where \(S\) is smooth over the mainlobe, \(X(f)\) is circular complex Gaussian, so \(\lvert X\rvert^2\) is exponential:

\[ P(f) \approx S(f)\, \frac{\chi^2_2}{2}, \qquad \mathrm{E}\{P\} = S, \qquad \mathrm{var}\{P\} = S^2, \qquad \frac{\mathrm{std}}{\mathrm{mean}} = 1 \]
(27.4)
\(P(f)\)
one periodogram bin
\(S(f)\)
true spectrum at that frequency
\(\chi^2_2\)
chi-squared variable with 2 degrees of freedom (real and imaginary parts); \(\chi^2_2/2\) is a unit-mean exponential
\(\mathrm{E},\ \mathrm{var},\ \mathrm{std}\)
mean, variance and standard deviation over realizations
\(\approx\)
valid away from the band edges and for a smooth \(S\)
Bins spaced \(f_s/N\) are approximately independent (rectangular window). A longer record gives more bins, not less noise per bin: the periodogram is not a consistent estimator.
\[ 10 \log_{10}(P/S):\quad \text{mean} = -\frac{10\,\gamma_E}{\ln 10} = -2.51\ \mathrm{dB}, \qquad \text{std} = \frac{10}{\ln 10}\cdot\frac{\pi}{\sqrt{6}} = 5.57\ \mathrm{dB} \]
(27.5)
\(10\log_{10}(P/S)\)
periodogram bin relative to the true spectrum, dB
\(P,\ S\)
as in (27.4)
\(\gamma_E\)
Euler constant, 0.5772
\(\ln\)
natural logarithm
\(\pi\)
3.14159
\(\gamma_E = 0.5772\) (Euler). The same 5.57 dB as log-compressed speckle (Module 08, Module 26): a periodogram bin and a speckle intensity pixel are both exponential random variables.
Demo 27.1

Periodogram of speckle RF: more samples, same scatter

SimulationWhite Gaussian noise through a Gaussian transfer function (FFT filtering), periodograms of 30–120 records, JS.
one periodogrammean over recordstrue S(f)e−x (χ²2/2)
Try: step \(N\) from 32 to 4096. The bin width shrinks from 1.25 MHz to 0.01 MHz, but \(\mathrm{std}(P/S)\) stays at 1 and the dB scatter at 5.6 dB. Switch to rectangular and look below −40 dB: the leakage floor is the window sidelobes convolved with the steep Gaussian band edge (bias); Hann pushes it down. Process: white Gaussian noise filtered by a 5 MHz, 60 % (−6 dB, amplitude) Gaussian transfer function, \(f_s = 40\) MHz.
Ultrasound link

The spectrum of one RF gate is itself speckle, in frequency. A 4 mm gate at 5 MHz gives a spectrum whose every bin has 100 % standard deviation. QUS gets its averaging from \(K\) laterally independent A-lines (separated by at least a beam width) and from the bandwidth, never from the gate alone. When a paper reports a spectrum, ask how many independent looks went into it.

Pitfall

Zero padding (Module 04) interpolates the same noisy spectrum onto a finer grid; smoothing adjacent padded bins averages correlated values and reduces variance far less than the number of bins suggests. Count independent bins at spacing ≈ ENBW, not FFT bins.

27.2Averaging: Bartlett, Welch and honest error bars

Bartlett. Split \(N\) samples into \(K\) non-overlapping segments of length \(M = N/K\) and average their periodograms. Variance falls to \(S^2/K\); resolution worsens to \(f_s/M = K f_s/N\). For fixed \(N\) the product is fixed:

\[ \frac{\mathrm{var}\{\hat S\}}{S^2} \times \Delta f_{\mathrm{res}} \approx \frac{1}{K} \times \frac{K f_s}{N} = \frac{f_s}{N} = \frac{1}{T} \]
(27.6)
\(K\)
number of averaged segments (looks)
\(\Delta f_{\mathrm{res}}\)
resolution cell, MHz
\(T\)
record duration \(N/f_s\), µs
\(\hat S\)
averaged spectrum estimate, as in (27.8)
\(S\)
true spectrum
\(\mathrm{var}\)
variance
\(N,\ f_s\)
record length (samples) and sampling rate (MHz), as in (27.1)
\(\times\)
multiplication (variance-resolution product)
Variance × resolution × record duration ≈ constant (up to a window-dependent factor of order 1). Only more data, or a smoother true spectrum, beats it.

Welch. Use windowed segments that overlap by \(M - D\) samples (step \(D\)). Tapering discards data at the segment edges; overlap recovers it. For a Gaussian process with \(S\) smooth over the window mainlobe, Welch (1967) gives

\[ \frac{\mathrm{var}\{\hat S\}}{S^2} = \frac{1}{N_b}\left[ 1 + 2 \sum_{m=1}^{N_b-1} \left(1 - \frac{m}{N_b}\right) \rho(mD)^2 \right] \equiv \frac{1}{K_{\mathrm{eff}}}, \qquad \rho(s) = \frac{\sum_n w[n]\, w[n+s]}{\sum_n w[n]^2} \]
(27.7)
\(N_b\)
number of segments
\(M, D\)
segment length and step, samples
\(\rho(s)\)
window overlap correlation at shift \(s\)
\(K_{\mathrm{eff}}\)
effective number of independent looks
\(\hat S,\ S\)
as in (27.6)
\(m\)
segment shift index
\(w[n],\ n\)
window and sample index, as in (27.2)
\(\sum_n\)
sum over window samples
\(\equiv\)
defined as
Hann at 50 % overlap: \(\rho(M/2) = 1/6\), so \(\rho^2 = 0.028\) and the segments are almost independent. \(N = 1024\), \(M = 128\): \(N_b = 15\), \(K_{\mathrm{eff}} = 14.3\) (m27_periodogram_welch.m measures var = 0.074 vs 0.070 predicted).

Confidence intervals. \(\hat S/S\) is approximately \(\chi^2_\nu/\nu\) with \(\nu = 2K_{\mathrm{eff}}\) equivalent degrees of freedom, so a 95 % interval for \(S\) is

\[ \left[\, \frac{\nu \hat S}{\chi^2_{\nu,\,0.975}}\,,\ \ \frac{\nu \hat S}{\chi^2_{\nu,\,0.025}} \,\right] \]
(27.8)
\(\nu\)
equivalent degrees of freedom, \(2K_{\mathrm{eff}}\)
\(\chi^2_{\nu,p}\)
\(p\)-quantile of the \(\chi^2_\nu\) distribution
\(\hat S\)
averaged spectrum estimate
\(0.975,\ 0.025\)
quantile levels for a 95 % confidence interval
\([\cdot,\ \cdot]\)
lower and upper 95 % limits of the true \(S\)
\(\nu = 2\) (one raw periodogram): \(S\) lies between \(\hat S\) − 5.67 dB and \(\hat S\) + 15.97 dB, a 21.6 dB span. \(\nu = 28.5\) (the Welch case above): −1.99 to +2.60 dB, a 4.6 dB span. The interval is constant in dB, which is why spectra with error bars are plotted on a log scale. Empirical coverage in the MATLAB script: 0.947 and 0.941.
Demo 27.2

Welch: variance against resolution on a rippled spectrum

SimulationGaussian RF with a two-echo transfer function (FFT filtering), Welch averages of 30 records, JS; the 95 % band uses the Wilson–Hilferty χ² quantile approximation.
one Welch estimate95 % χ² bandmean of 30 (bias)true S(f)
Try: \(N = 1024\) samples (25.6 µs). The process has two echoes \(\tau = 2\) µs apart (amplitude ratio 0.5), which writes a ripple of period \(1/\tau = 0.5\) MHz and 9.5 dB depth onto the spectrum. Shrink \(M\): the band narrows (variance \(\propto 1/K_{\mathrm{eff}}\)) and the mean estimate flattens. The ripple survives only for \(M\) well above \(f_s\tau = 80\) samples: at \(M = 128\) the Hann ENBW (0.47 MHz) already matches the ripple period and the mean is nearly flat; at \(M = 64\) a segment cannot even contain both echoes. That is bias, not noise. Compare 0 % and 50 % overlap with Hann: almost twice the segments for nearly the same independence. The measured variance exceeds the Welch formula here because the echo pair correlates the data across segments; the formula assumes \(S\) flat over the window.
Intuition

Two echoes separated by \(\tau\) multiply the spectrum by \(\lvert 1 + a\, e^{-j2\pi f\tau}\rvert^2\), a cosine ripple with period \(1/\tau\). This is the spectral signature of quasi-periodic scatterer spacing (mean scatterer spacing estimation looks for exactly this ripple). The estimator can only see it if its segment, or gate, is longer than \(\tau\), which ties spectral resolution directly to the axial extent of the gate.

27.3Multitaper spectra: Slepian sequences

Welch throws away the gate edges and uses only part of the record per segment. Thomson's multitaper method uses the whole record \(K\) times, each time with a different orthogonal taper chosen to leak as little as possible outside a band \(\pm W\). The best taper solves a concentration problem (Slepian, 1978):

\[ \lambda = \frac{\int_{-W}^{W} \lvert V(f)\rvert^2\, df}{\int_{-1/2}^{1/2} \lvert V(f)\rvert^2\, df} \to \max \quad\Longleftrightarrow\quad \mathbf{A}\mathbf{v} = \lambda \mathbf{v}, \qquad A_{mn} = \frac{\sin\!\left(2\pi W(m-n)\right)}{\pi(m-n)} \]
(27.9)
\(\mathbf{v}\)
taper of length \(N\); \(V(f)\) its DTFT
\(W\)
half-bandwidth, cycles per sample
\(\lambda\)
fraction of taper energy inside \(\pm W\) (concentration)
\(\mathbf{A}\)
\(N \times N\) sinc kernel matrix
\(V(f)\)
DTFT of the taper, as above
\(f\)
frequency in cycles per sample
\(\int_{-W}^{W}\)
energy inside the band \(\pm W\)
\(\int_{-1/2}^{1/2}\)
total energy (Nyquist band)
\(A_{mn}\)
element (m, n) of \(\mathbf{A}\); \(m, n = 0 \dots N-1\)
\(\Longleftrightarrow\)
equivalent to
\(\sin\)
sine; the kernel is the sinc function sampled at integer lags
\(f\) and \(W\) in cycles per sample. The eigenvectors \(\mathbf{v}_k\) are the discrete prolate spheroidal sequences (DPSS), orthonormal, ordered by concentration \(\lambda_0 > \lambda_1 > \dots\). About \(2NW\) eigenvalues are close to 1 (the Shannon number), then they fall off a cliff.
\[ T_{nn} = \left(\frac{N-1}{2} - n\right)^2 \cos 2\pi W, \qquad T_{n,n+1} = \frac{(n+1)(N-n-1)}{2} \]
(27.10)
\(T_{nn}\)
diagonal entries of the tridiagonal matrix \(\mathbf{T}\)
\(T_{n,n+1}\)
off-diagonal entries of \(\mathbf{T}\)
\(n\)
row index, \(n = 0 \dots N-1\)
\(N,\ W\)
taper length (samples) and half-bandwidth (cycles/sample), as in (27.9)
\(\mathbf{T}\)
matrix that commutes with \(\mathbf{A}\) of (27.9), so it has the same eigenvectors (the DPSS tapers)
This tridiagonal matrix commutes with \(\mathbf{A}\), so it has the same eigenvectors and is far better conditioned numerically (\(\mathbf{A}\)'s eigenvalues crowd near 1). Its eigenvectors with the largest eigenvalues are the DPSS. m27_multitaper_dpss.m builds them with eig and checks them against \(\mathbf{A}\).
\[ \hat S_{\mathrm{MT}}(f) = \frac{1}{K} \sum_{k=0}^{K-1} \left\lvert \sum_n v_k[n]\, x[n]\, e^{-j2\pi f n} \right\rvert^2, \qquad K \approx 2NW - 1, \qquad \mathrm{var} \approx \frac{S^2}{K}, \qquad \text{resolution } 2W \]
(27.11)
\(v_k[n]\)
\(k\)-th DPSS taper
\(K\)
number of tapers (eigenspectra averaged)
\(NW\)
time-bandwidth product, dimensionless
\(\hat S_{\mathrm{MT}}(f)\)
multitaper spectrum estimate
\(x[n],\ f\)
as in (27.1); here f in cycles per sample
\(N,\ W\)
as in (27.9)
\(\mathrm{var}\)
variance of the estimate
\(S\)
true spectrum
\(2W\)
resolution (full bandwidth of the spectral window), cycles per sample
Each eigenspectrum is approximately \(\chi^2_2\) and they are nearly uncorrelated for a locally flat spectrum, so \(\nu \approx 2K\). \(N = 128\) (2.5 mm gate at 40 MHz), \(NW = 4\): \(2W = 2NW f_s/N = 2.5\) MHz, \(K = 7\), \(\lambda_6 = 0.937\), \(\lambda_7 = 0.699\). Thomson's adaptive weighting downweights the leaky high-order eigenspectra where the spectrum is steep.
Demo 27.3

DPSS tapers, their spectral windows, and the multitaper estimate of a short gate

SimulationDPSS from the tridiagonal eigenproblem (27.10) by implicit QL, multitaper spectra of 60 simulated Gaussian gates, JS.
v0 / multitapervK−1 / true Srectangular periodogram
Try: at \(NW = 4\) raise \(K\) from 7 to 8: \(\lambda_{K-1}\) drops from 0.94 to 0.70 and the out-of-band leakage of the last taper rises to −5 dB (at \(K = 9\) the ninth taper has \(\lambda_8 = 0.30\)). Raise \(NW\): variance falls as \(1/K\) but the \(2W\) smoothing biases the peak of this short-gate spectrum low (mean \(\hat S/S\) below 1). Look above 10 MHz: the rectangular periodogram flattens near −40 dB while the true spectrum keeps falling. With \(K = 2NW - 1\) the unweighted multitaper average leaks too, because its last tapers (\(k\) = 4–6, \(\lambda\) = 0.999–0.94) put \(1 - \lambda\) = −32 to −12 dB of their energy outside \(\pm W\); drop \(K\) to 3 and the floor falls by tens of dB. Thomson's adaptive weights make this choice per frequency. Variance is measured over 60 independent gates.
Pitfall

Multitaper is not a free lunch for a smooth speckle spectrum: for the same spectral-window bandwidth, Hann–Welch and multitaper give nearly the same variance (m27_multitaper_dpss.m: variance × bandwidth × \(T\) = 1.36 vs 1.39). Its real advantages are leakage control over a large dynamic range (with low-order tapers or adaptive weighting), use of the full gate, and a principled bandwidth knob. Every extra high-order taper buys variance with leakage, and beyond \(K = 2NW - 1\) the price is steep.

27.4Parametric (AR) spectra for short gates

A Fourier estimator cannot resolve finer than about \(1/T\). An autoregressive model assumes the data are white noise through an all-pole filter and estimates the few pole parameters instead of the spectrum bin by bin:

\[ x[n] = -\sum_{k=1}^{p} a_k\, x[n-k] + e[n], \qquad S_{\mathrm{AR}}(f) = \frac{\sigma_e^2}{\left\lvert 1 + \sum_{k=1}^{p} a_k\, e^{-j2\pi f k/f_s} \right\rvert^2} \]
(27.12)
\(p\)
model order
\(a_k\)
AR (prediction) coefficients
\(e[n]\)
white driving noise of variance \(\sigma_e^2\)
\(S_{\mathrm{AR}}(f)\)
AR power spectrum
\(\sigma_e^2\)
variance of \(e[n]\)
\(f,\ f_s\)
frequency and sampling rate, MHz
\(x[n],\ n\)
signal and sample index, as in (27.1)
\(k\)
coefficient index, 1 … p
\(j\)
imaginary unit
\[ \text{Yule–Walker:}\quad \mathbf{R}\,\mathbf{a} = -\mathbf{r}, \qquad R_{ij} = \hat r\!\left[\lvert i-j\rvert\right], \qquad \text{solved by Levinson–Durbin in } O(p^2) \]
(27.13)
\(\mathbf{R}\)
\(p \times p\) Toeplitz matrix of estimated autocorrelations
\(\mathbf{a}\)
vector of AR coefficients \((a_1 \dots a_p)\), as in (27.12)
\(\mathbf{r}\)
vector \((\hat r[1] \dots \hat r[p])\)
\(\hat r[l]\)
estimated autocorrelation at lag \(l\), as in (27.1)
\(i,\ j\)
row and column indices
\(O(p^2)\)
operation count of the Levinson–Durbin recursion
With the biased estimate \(\hat r[l] = (1/N)\sum x[n]\,x[n+l]\), which implicitly windows the data with a triangle and smears close peaks. The model is always stable.
\[ \text{Burg:}\quad k_m = \frac{-2 \sum_n f_{m-1}[n]\, b_{m-1}[n-1]}{\sum_n \left( f_{m-1}[n]^2 + b_{m-1}[n-1]^2 \right)}, \qquad \sigma_m^2 = \sigma_{m-1}^2 \left(1 - k_m^2\right) \]
(27.14)
\(f_m, b_m\)
forward and backward prediction errors at order \(m\)
\(k_m\)
reflection coefficient, \(\lvert k_m\rvert \le 1\)
\(\sigma_m^2\)
prediction-error power at order \(m\)
\(n\)
sample index
\(m\)
recursion order
The coefficients follow the Levinson update \(\mathbf{a}^{(m)} = [\mathbf{a}^{(m-1)}; 0] + k_m\, \mathrm{flip}([\mathbf{a}^{(m-1)}; 0])\). \(\lvert k_m\rvert \le 1\) guarantees stability, and no window is assumed outside the data.
\[ \mathrm{AIC}(p) = N \ln \hat\sigma_p^2 + 2p \]
(27.15)
\(\mathrm{AIC}(p)\)
Akaike information criterion for order p; the minimum selects the order
\(N\)
number of data samples
\(\hat\sigma_p^2\)
estimated prediction-error power of the order-p fit, as in (27.14)
\(p\)
model order, as in (27.12)
\(\ln\)
natural logarithm
Akaike (1974). Picks a reasonable order for the two-tone case (median \(p = 16\) for \(N = 64\) at 20 dB SNR); it tends to over-fit for long records, where MDL (penalty \(p \ln N\)) is more conservative.

Numbers from m27_ar_burg.m, two tones 0.4 MHz apart in a 64-sample gate (Fourier resolution \(f_s/N = 0.625\) MHz, 20 dB SNR per tone): probability of two correctly placed peaks 0.37 / 0.66 / 0.75 for Burg at \(p\) = 12 / 16 / 24, 0.00 / 0.09 / 0.29 for Yule–Walker, 0.02 for the Hann periodogram.

Demo 27.4

Burg vs Yule–Walker vs periodogram on a short gate

SimulationTwo tones in white Gaussian noise; Burg, Levinson–Durbin and Hann-DTFT spectra of each gate; 100 gates per readout, JS.
Hann periodogramBurg AR(p)Yule–Walker AR(p)true tones / AIC
Try: at \(N = 64\) the tones are inside one Fourier resolution cell. Raise \(p\): Burg splits them first, Yule–Walker much later. Push \(p\) near \(N/2\): spurious peaks appear (P(>2 peaks) rises). Drop SNR to 5 dB: AR resolution collapses faster than the periodogram's, because the poles retreat from the unit circle. The probabilities come from 100 independent gates; peaks count if they rise within 20 dB of the maximum.
Intuition

For a smooth speckle spectrum, what matters is low-order structure, so AR helps the spectral moments more than the resolution. In m27_ar_burg.m the centroid of a 64-sample speckle gate has std 0.36 MHz from a Hann periodogram and 0.29 MHz from Burg at any order from 4 to 16. The reason: the one-sided centroid equals \(f_s/4 - (2f_s/\pi^2) \sum_{\text{odd } l} \rho[l]/l^2\), dominated by lag 1, and every Burg model of order \(\ge l\) reproduces the same \(\rho[0..l]\).

Pitfall

AR peak heights are not powers (the area under the peak is), Burg peaks shift and split with initial phase at high SNR, and an AR spectrum has no \(\chi^2\) error bar. Use AR for locating features and for moments on very short gates; use averaged Fourier spectra when you need calibrated power, as in BSC estimation.

27.5From spectrum to tissue: the QUS signal model and the reference phantom

For incoherent scattering from many sub-resolution scatterers, the expected power spectrum of a gate centered at depth \(z\) factorizes into system, beam, propagation and tissue terms:

\[ S(f, z) = G(f)\cdot D(f, z)\cdot \exp\!\left( -4 \int_{0}^{z} \alpha(f, z')\, dz' \right) \cdot \mathrm{BSC}(f)\cdot C_{\mathrm{gate}} \]
(27.16)
\(G(f)\)
transmit–receive electromechanical response (power)
\(D(f,z)\)
diffraction and focusing, including the frequency-dependent beam cross-section
\(\alpha(f,z)\)
amplitude attenuation coefficient, Np/cm
\(\mathrm{BSC}(f)\)
backscatter coefficient, sr⁻¹ cm⁻¹
\(C_{\mathrm{gate}}\)
gate length, window energy, gain and TGC
\(S(f,z)\)
measured power spectrum of the gate centred at depth z
\(f,\ z\)
frequency (MHz) and depth (cm)
\(z^{\prime}\)
integration variable (depth), cm
\(\exp(\cdot)\)
exponential function; the factor 4 is two-way travel (2) times power (2)
The exponential is the two-way power loss: the 4 = 2 for the round trip × 2 for power.

Reference phantom method (Yao, Zagzebski & Madsen, 1990). Scan a phantom with known \(\alpha_r(f)\) and \(\mathrm{BSC}_r(f)\) with the same transducer, focus, gain, TGC and gate. \(G\), \(D\) and \(C_{\mathrm{gate}}\) divide out:

\[ \frac{S_s(f, z)}{S_r(f, z)} = \frac{\mathrm{BSC}_s(f)}{\mathrm{BSC}_r(f)}\, \exp\!\left( -4\, [\alpha_s(f) - \alpha_r(f)]\, z \right) \qquad (\text{homogeneous media}) \]
(27.17)
\(S_s, S_r\)
sample and reference gate spectra at the same depth and settings
\(\alpha_s, \alpha_r\)
sample and reference attenuation, Np/cm
\(\mathrm{BSC}_s,\ \mathrm{BSC}_r\)
sample and reference backscatter coefficients, sr⁻¹ cm⁻¹
\(f,\ z\)
frequency (MHz) and depth (cm)
This one equation gives both QUS families: its depth dependence gives \(\alpha_s\); its frequency shape after attenuation compensation gives \(\mathrm{BSC}_s\). It assumes the same sound speed in sample and reference (otherwise the focus and \(D\) differ), the same overlying attenuation for BSC (or an explicit correction), and incoherent scattering (no specular interfaces in the gate).
Schematic Reference phantom method: what must match, what divides outSame probe, focus, gain, TGC and gate in both scans; the ratio at equal depth keeps only the BSC ratio and the attenuation difference. Typical numbers; not to scale.
Sample (tissue)Reference phantom (known αᵣ, BSCᵣ) linear probe, 5 MHz, F/2 same probe, focus, gain, TGC focus 2 cm depth z gate 4 mm (Hann)same gate, same z αₛ(f), BSCₛ(f), c ≈ 1540 m/sαᵣ(f) = 0.5 dB/cm/MHz, BSCᵣ(f) ∝ f⁴, same c Sₛ/Sᵣ = [G·D·C·BSCₛ·e^(−4αₛz)] / [G·D·C·BSCᵣ·e^(−4αᵣz)] = (BSCₛ/BSCᵣ)·e^(−4(αₛ−αᵣ)z) Breaks if c differs (focus moves, D no longer cancels), if overlying layers differ (BSC),or if αᵣ, BSCᵣ are wrong (they pass straight into the result).

Spectral moments. Centroid \(f_c = \int f S\, df / \int S\, df\) and spectral variance \(\sigma^2 = \int (f - f_c)^2 S\, df / \int S\, df\). For linear attenuation \(\alpha(f) = \alpha_0 f\), differentiate \(f_c\) with respect to depth:

\[ \frac{\partial f_c}{\partial z} = -4\, \alpha_0\, \sigma^2(z) \qquad (\text{exact, any spectral shape}) \]
(27.18)
\(f_c\)
power-spectrum centroid, MHz
\(\sigma^2(z)\)
power-spectrum variance at depth \(z\), MHz²
\(\alpha_0\)
attenuation slope, Np/cm/MHz
\(z\)
depth, cm
\(\partial f_c/\partial z\)
rate of change of the centroid with depth, MHz/cm
Because \(\partial S/\partial z = -4\alpha_0 f S\). For a Gaussian power spectrum \(\sigma^2 = \sigma_p^2\) is constant, so the downshift is linear: \(f_c = f_0 - 4\alpha_0\sigma_p^2 z\), the same as Module 08's \(-2\alpha_0\sigma_f^2 z\) with \(\sigma_f = \sqrt{2}\,\sigma_p\) the amplitude-spectrum std. For other shapes \(\sigma^2\) changes with depth and the downshift curves (card 4 below).
Ultrasound link

Commercial attenuation-imaging and controlled-attenuation-parameter style tools report \(\alpha\) in dB/cm/MHz (or dB/m at a fixed frequency) and must solve the same normalization problem; some use a built-in reference (a stored phantom or system model) rather than a phantom scan per exam. When comparing values across systems, check the frequency band, the depth range of the ROI and whether diffraction was compensated: in the toy model of m27_attenuation_sld.m an uncompensated focus moves \(\alpha\) from 0.70 to 0.96 dB/cm/MHz.

27.6Attenuation from spectra: log difference, regression and centroid

Write the log spectral ratio of sample to reference, in nepers of power:

\[ \begin{aligned} y(f; z) &= \ln S_s(f, z) - \ln S_r(f, z) \\ &= \ln\frac{\mathrm{BSC}_s}{\mathrm{BSC}_r}(f) - 4\, \Delta\alpha(f)\, z + \varepsilon(f, z), \qquad \Delta\alpha = \alpha_s - \alpha_r \ \ (\mathrm{Np/cm}) \end{aligned} \]
(27.19)
\(y(f;z)\)
log spectral ratio, Np of power
\(\Delta\alpha(f)\)
attenuation difference, Np/cm
\(\varepsilon(f,z)\)
log-spectral noise (mean \(\psi(K) - \ln K\), variance \(\psi_1(K)\))
\(S_s,\ S_r\)
as in (27.17)
\(\mathrm{BSC}_s,\ \mathrm{BSC}_r\)
as in (27.17)
\(\ln\)
natural logarithm
\(\alpha_s,\ \alpha_r\)
as in (27.17)
\(f,\ z\)
frequency (MHz) and depth (cm)

Spectral log difference (SLD), two gates at \(z_1 < z_2\), \(\Delta z = z_2 - z_1\). The BSC term cancels if the medium is homogeneous between them:

\[ \begin{aligned} y(f; z_1) - y(f; z_2) &= 4\, \Delta\alpha(f)\, \Delta z \\ \Rightarrow\ \ \text{with } \Delta\alpha = \Delta\alpha_0 f:\quad \hat\alpha_0 &= \alpha_{r,0} + \frac{\text{slope vs } f}{4\, \Delta z}\ \ [\mathrm{Np/cm/MHz}], \qquad \times 8.686 \to \mathrm{dB/cm/MHz} \end{aligned} \]
(27.20)
\(z_1, z_2\)
proximal and distal gate centers, cm; \(\Delta z = z_2 - z_1\)
\(\alpha_{r,0}\)
reference attenuation slope, Np/cm/MHz
\(\hat\alpha_0\)
estimated sample attenuation slope, Np/cm/MHz
8.686
dB per Np (\(20/\ln 10\))
\(y(f;z_1),\ y(f;z_2)\)
log spectral ratio at the two gates, as in (27.19)
\(\Delta z\)
gate separation \(z_2 - z_1\), cm
\(\Delta\alpha_0\)
attenuation-slope difference, Np/cm/MHz (so \(\Delta\alpha = \Delta\alpha_0 f\))
\(\text{slope vs } f\)
slope of the fitted line of the log-ratio difference versus frequency, Np/MHz
\(\times 8.686\)
conversion from Np to dB

Reference phantom regression (spectral difference): at each frequency, regress \(y(f; z)\) on \(z\) over all gates in the ROI to get \(-4\Delta\alpha(f)\), then fit \(\Delta\alpha(f) = \Delta\alpha_0 f\). Centroid downshift: regress \(f_c(z)\) on \(z\) and divide by \(-4\sigma_p^2\); subtract the reference phantom's centroid slope to remove the system and diffraction terms to first order.

\[ \begin{aligned} \ln(\text{power ratio}) &= -4\, \alpha_{\mathrm{Np}}\, f\, \Delta z \\ \Longleftrightarrow\quad 10 \log_{10}(\text{power ratio}) &= -2\, \alpha_{\mathrm{dB}}\, f\, \Delta z\ \ \mathrm{dB} \\ \Longleftrightarrow\quad \ln(\text{amplitude ratio}) &= -2\, \alpha_{\mathrm{Np}}\, f\, \Delta z \end{aligned} \]
(27.21)
\(\alpha_{\mathrm{Np}}\)
attenuation slope, Np/cm/MHz
\(\alpha_{\mathrm{dB}}\)
the same slope in dB/cm/MHz, \(8.686\,\alpha_{\mathrm{Np}}\)
\(f\)
frequency, MHz
\(\Delta z\)
path length in the medium, cm
\(\log_{10},\ \ln\)
base-10 and natural logarithms
1 Np = 8.686 dB (\(20/\ln 10\)). 0.5 dB/cm/MHz = 0.0576 Np/cm/MHz. At 5 MHz over \(\Delta z = 2\) cm: ln power ratio −2.30, which is −10.0 dB of power (2 × 0.5 × 5 × 2). Use amplitude spectra and the factor is 2, not 4; use dB and the 8.686 is already inside \(\alpha_{\mathrm{dB}}\). Most published factor-of-2 errors come from mixing these.

Variance. Averaging \(K\) independent lines makes \(\ln\hat S - \ln S\) a log of a mean of \(K\) exponentials: mean \(\psi(K) - \ln K\) (cancels in a difference when both gates use the same \(K\)), variance \(\psi_1(K) \approx 1/K + 1/(2K^2)\) (trigamma). The band \(B\) holds about \(BT/\gamma_w\) independent frequency samples (\(T\) = gate duration). A least-squares slope over \(B\) then gives, for the two-gate SLD,

\[ \sigma(\hat\alpha_0) \approx \frac{1}{\Delta z} \sqrt{ \frac{1.5\, \gamma_w\, \psi_1(K)}{B^3\, T} } \qquad [\mathrm{Np/cm/MHz};\ B \text{ in MHz},\ T \text{ in } \mu\mathrm{s},\ \Delta z \text{ in cm}] \]
(27.22)
\(\gamma_w\)
window bin-correlation factor: 0.93 (rectangular gate), 1.75 (Hann)
\(\psi_1(K)\)
trigamma function, variance of the log of a \(K\)-look average
\(B\)
analysis bandwidth, MHz
\(T\)
gate duration, µs
\(\sigma(\hat\alpha_0)\)
standard deviation of the slope estimate, Np/cm/MHz
\(\hat\alpha_0\)
estimated slope, as in (27.20)
\(\Delta z\)
gate separation, cm
\(K\)
number of independent spectral looks per gate
\(\gamma_w\) comes from Monte Carlo on white noise; it measures the bin correlation of the windowed spectrum. \(K = 32\), Hann 4 mm gates (\(T = 5.2\) µs), \(B = 4\) MHz, \(\Delta z = 2\) cm: 0.069 dB/cm/MHz predicted, 0.067 measured over 128 ROIs in m27_attenuation_sld.m.
Demo 27.5

Attenuation by the spectral log difference: bias, spread and the factor 4

Approx. modelMultiplicative speckle model: gate spectrum = S(f, z) × periodogram of a Hann-windowed white Gaussian gate averaged over K lines (exact when S is smooth over the 2/T window); noise-free reference; no electronic noise or diffraction; 300 ROIs, JS.
y(f; z₁) − y(f; z₂), one ROI, and its fittruthα̂ over 300 ROIs
Try: halve \(\Delta z\) or the band: \(\sigma\) doubles and grows 2.8× respectively (\(1/\Delta z\), \(B^{-3/2}\)). Quadruple \(K\): \(\sigma\) halves. Pick the "error" conversion: the estimate of \(\alpha_s - \alpha_r\) doubles (\(\alpha_r\) = 0.5 dB/cm/MHz). Model: gate spectra are \(S(f, z)\) times the periodogram of a Hann-windowed white Gaussian gate averaged over \(K\) lines (the multiplicative speckle model, exact when \(S\) is smooth over the \(2/T\) window); the reference is averaged over many lines and treated as noise-free; \(z_1\) = 2 cm.
Demo 27.6

Centroid downshift vs log-ratio regression on the same data

Approx. modelToy focal gain D = 1/(1 + ((z − z_F)/z_R)²), z_R = 3.5λF#², times the multiplicative speckle model; Gaussian or two-lobe transducer spectrum; BSC ∝ fⁿ; 150 ROIs, JS.
sample centroid (one ROI, expected)reference centroidnominal −4α₀σp²ztrue α
Try: with focusing on, the raw centroid reads about 0.87 for a true 0.70 (the focal gain is frequency dependent, so it shifts the centroid with depth); subtracting the reference centroid slope brings it to about 0.67. Turn focusing off: with \(n = 4\) the raw centroid reads about 0.60, because \(f^4\) narrows the spectral variance from the nominal 0.81 to about 0.71 MHz², exactly the ratio \(\partial f_c/\partial z = -4\alpha_0\sigma^2\) predicts. Set \(n = 0\): the raw centroid is right, and subtracting a reference whose BSC goes as \(f^4\) now makes it worse (about 0.75), because the subtraction is exact only when both spectra have the same variance. Switch to the asymmetric spectrum: both centroid estimators fail (0.4 to 0.6), while the log-ratio regression, which assumes nothing about spectral shape, stays at 0.70. ROI: 8 contiguous 4 mm Hann gates, 1.7–4.5 cm; focus 2 cm, F/2; reference \(\alpha_r\) = 0.5 dB/cm/MHz, \(\mathrm{BSC}_r \propto f^4\); 150 ROIs.
Bias sources, ranked roughly by size in practice

(1) Diffraction and sound-speed mismatch between sample and reference. (2) Electronic noise: it flattens deep, high-frequency spectra and biases \(\alpha\) low; restrict the band to where SNR exceeds ~10 dB or subtract the noise power spectrum. (3) BSC that changes with depth inside the ROI (layered or heterogeneous tissue): the log-ratio method reads it as attenuation. (4) Specular reflectors and coherent components in the gate. (5) Unequal effective number of looks at the two depths, which breaks the cancellation of \(\psi(K) - \ln K\). Window smoothing and intra-gate attenuation are second order for gates of a few millimetres.

27.7Backscatter coefficient, form factors, ESD and Lizzi parameters

The BSC is the differential scattering cross-section per unit volume at 180°, in sr⁻¹ cm⁻¹ (or m⁻¹ sr⁻¹). From the reference ratio, after compensating the attenuation difference over the full path to depth \(z\):

\[ \mathrm{BSC}_s(f) = \mathrm{BSC}_r(f)\, \frac{S_s(f, z)}{S_r(f, z)}\, \exp\!\left( 4\, [\alpha_s(f) - \alpha_r(f)]\, z \right) \]
(27.23)
\(\mathrm{BSC}_s,\ \mathrm{BSC}_r\)
sample and reference backscatter coefficients, sr⁻¹ cm⁻¹
\(S_s,\ S_r\)
sample and reference gate spectra, as in (27.17)
\(\alpha_s,\ \alpha_r\)
sample and reference attenuation coefficients, Np/cm
\(f,\ z\)
frequency (MHz) and depth (cm)
\(\exp(\cdot)\)
exponential function; compensates the differential two-way power attenuation

For weak scatterers small compared with \(\lambda\) the BSC rises as \(f^4\) (Rayleigh). Finite size multiplies it by a form factor \(F\), the normalized spatial autocorrelation of the impedance fluctuations in the wavenumber domain. The Gaussian form factor (Insana, Wagner, Brown & Hall, 1990) is the workhorse:

\[ \mathrm{BSC}(f) \propto k^4 \exp\!\left(-0.827\, k^2 a^2\right), \qquad k = \frac{2\pi f}{c}, \qquad \mathrm{ESD} = 2a \]
(27.24)
\(k\)
wavenumber, rad/µm (with \(c\) = 1540 µm/µs and \(f\) in MHz)
\(a\)
effective scatterer radius, µm
\(\mathrm{ESD}\)
effective scatterer diameter, µm
0.827
constant of the Gaussian form factor \(F = \exp(-0.827\,k^2a^2)\)
\(\mathrm{BSC}(f)\)
backscatter coefficient, as in (27.16)
\(f\)
frequency, MHz
\(c\)
sound speed, µm/µs
\(\propto\)
proportional to
Taking logs, \(\ln[\mathrm{BSC}_s/\mathrm{BSC}_r]\) (with a Rayleigh reference, \(\mathrm{BSC}_r \propto f^4\)) is linear in \(k^2\) with slope \(-0.827a^2\): an ordinary least-squares problem. The intercept is proportional to the effective acoustic concentration (number density × squared impedance contrast).
\[ n_{\mathrm{local}}(ka) = \frac{d \ln \mathrm{BSC}}{d \ln f} = 4 - 1.654\, (ka)^2, \qquad k^4F \text{ peaks at } ka = \sqrt{2/0.827} = 1.555 \]
(27.25)
\(n_{\mathrm{local}}\)
local power-law exponent of the BSC versus frequency
\(ka\)
wavenumber times effective radius, dimensionless
\(\mathrm{BSC}\)
as in (27.24)
\(f\)
frequency, MHz
\(k^4F\)
Rayleigh \(k^4\) law times the form factor \(F = \exp(-0.827\,k^2a^2)\)
\(d\ln\mathrm{BSC}/d\ln f\)
logarithmic slope of the BSC with frequency
The spectral slope carries the size information only where \(n\) departs from 4. With a 3–7 MHz band, ESD = 20 µm spans \(ka\) = 0.12–0.29 (\(n \approx 3.9\) at 5 MHz): indistinguishable from Rayleigh. ESD = 100 µm spans \(ka\) = 0.61–1.43 (\(n = 2.3\) at 5 MHz): well identified.

Lizzi parameters (Lizzi et al., 1983). Fit a straight line to the calibrated, attenuation-compensated spectrum in dB over the band: the slope (dB/MHz) falls as scatterers grow; the intercept (dB at 0 MHz) depends on size, concentration and impedance contrast; the midband fit (the line's value at band center) combines them and is the most precise of the three because a least-squares line is best determined at the mean abscissa.

Demo 27.7

ESD from a Gaussian form factor: identifiability and attenuation leakage

Approx. modelGaussian form factor (27.24) times the multiplicative speckle model of Demo 27.5; 4 mm Hann gate at 3 cm; Rayleigh reference with equal attenuation; band assumed above the noise floor; 300 ROIs, JS.
ln(BSCs/BSCr), one ROI, and fittruthanalysis band in kaESD over 300 ROIs
Try: slide the ESD down to 20 µm at 3–7 MHz: the \(k^2\) slope is buried in noise, many fits return a negative \(a^2\) (no solution) and the survivors are biased high. Move to 8–15 MHz and the same scatterers become measurable. Add +0.10 dB/cm/MHz of compensation error at 3 cm: ESD 100 µm reads about 78 µm, because \(\exp(4\delta\alpha f z)\) is a linear-in-\(f\) tilt that the \(k^2\) fit partly absorbs. 4 mm Hann gate at \(z\) = 3 cm, Rayleigh reference with equal attenuation, same multiplicative speckle model as Demo 27.5 (the band is assumed to be above the noise floor).

27.8Gate, ROI and band: the precision budget

Attenuation imaging trades spatial resolution for precision more steeply than any other QUS parameter. Regress the log ratio over an ROI of axial length \(L\) tiled with \(n\) contiguous gates of length \(cT/2\). Each gate contributes a frequency slope with variance \(12\gamma_w\psi_1(K)/(B^3T)\); the depth regression over \(n = L/(cT/2)\) gates divides by \(\sum(z - \bar z)^2 \approx n^3(cT/2)^2/12\):

\[ \sigma(\hat\alpha_0) \approx \sqrt{ \frac{4.5\, \gamma_w\, \psi_1(K)\, c}{B^3\, L^3} } \qquad [\mathrm{Np/cm/MHz};\ c \text{ in cm}/\mu\mathrm{s}] \]
(27.26)
\(L\)
axial ROI length, cm
\(K\)
independent lines (lateral looks) in the ROI
\(B\)
analysis bandwidth, MHz
\(c\)
sound speed, 0.154 cm/µs
\(\gamma_w\)
1.75 for Hann gates, 0.93 for rectangular
\(\sigma(\hat\alpha_0)\)
standard deviation of the slope estimate, Np/cm/MHz
\(\hat\alpha_0\)
estimated slope, as in (27.20)
\(\psi_1(K)\)
trigamma function, as in (27.22)
\(\approx\)
valid when the ROI length \(L\) replaces the gate separation \(\Delta z\) of (27.22)
The gate duration \(T\) cancels: shorter gates give noisier spectra but more of them. Monte Carlo (2, 4, 8 mm gates, rectangular and Hann) agrees within 5 %. Hann, \(K = 32\), \(B = 4\) MHz, \(L = 3\) cm: 0.041 dB/cm/MHz. Overlapping Hann gates (75 %) recover the tapered edges and do somewhat better (0.030 dB/cm/MHz in m27_attenuation_sld.m for 1.5–4.5 cm). Two gates of length \(L/2\) at the ends give the constant 6 instead of 4.5, and two short gates far apart are worse still.

So precision scales as \(K^{-1/2} B^{-3/2} L^{-3/2}\): doubling the ROI length buys 2.8×, doubling the number of lines (the lateral ROI width) only 1.4×. The gate length is chosen for other reasons: long enough that the window bandwidth (≈ \(2/T\)) is small against spectral structure and that the log-average bias is the same at all depths; short enough that the medium is homogeneous inside it.

Schematic Gate and ROI layout: two-gate SLD versus regression over contiguous gatesOne attenuation ROI seen from the probe. Lines are spaced by about one beam width so that they count as independent looks. Typical numbers for 5 MHz, F/2; gates drawn wider than scale.
Two-gate SLD (27.20), (27.22)Regression over contiguous gates (27.26) gate at z₁gate at z₂ Δz = 2 cm cT/2 = 4 mmn = L/(cT/2)gates L = 3 cm (1.5–4.5 cm) K = 32 lines × ≈ 0.6 mm (λF# at 5 MHz, F/2) ≈ 2 cm same K lines; every gate's spectrum enters the depth fit

27.9Visual summary

Each card is a question you should answer in five seconds, the answer, and the picture to keep in mind. The slider rescales the gate in the first three cards: a longer gate always means a narrower spectral window and a finer ripple.

gate T = 1.00 → spectral window ∝ 1/T = 1.00

27.10How experts think about it

An expert sees a QUS estimate as a chain: a spectrum that is exponentially distributed per bin, averaged over a known number of independent looks, divided by a reference that removes everything except attenuation and backscatter, then fitted with a model whose parameters trade off against each other. Every number gets a degrees-of-freedom count, a unit check (Np or dB, amplitude or power, one-way or two-way) and a list of what was assumed to cancel.

A periodogram bin has 100 % standard deviation, however long the record.
std/mean = 1\(\chi^2_2/2\)5.57 dB in dB units

\(X(f)\) of a Gaussian process is circular complex Gaussian, so \(\lvert X\rvert^2\) is exponential with mean \(S\) and std \(S\) (27.4). A longer record makes more bins (spacing \(f_s/N\)), each just as noisy. m27_periodogram_welch.m: \(\mathrm{std}(P/S)\) = 1.00, 0.99, 1.00 for \(N\) = 64, 256, 1024.

SimulationDots: FFT periodograms of one white Gaussian record per N, JS. Lines: 1 and \(1/\sqrt{K}\) (closed form).
Near 0 and \(f_s/2\) the bin is \(\chi^2_1\) (variance \(2S^2\)); non-Gaussian data (few scatterers, specular echoes) and strong leakage change the distribution.
Averaging K looks buys variance 1/K and costs resolution K/N: variance × resolution × duration is fixed.
var \(\propto 1/K\)\(\Delta f \propto K/N\)Hann 50 %: \(\rho^2 = 0.028\)

Bartlett: \(K\) segments of \(N/K\) samples, var \(S^2/K\), resolution \(K f_s/N\) (27.6). Welch with Hann at 50 % overlap has segment correlation \(\rho = 1/6\), so it nearly doubles the segment count at almost no loss: \(N = 1024\), \(M = 128\) gives \(N_b = 15\) and \(K_{\mathrm{eff}} = 14.3\) (27.7).

Exact formula
Assumes \(S\) is smooth over the window mainlobe; with structure finer than \(f_s/M\) the "variance reduction" is really bias (the ripple in Demo 27.2 disappears).
Quote a spectrum with its degrees of freedom: the 95 % interval spans 21.6 dB at ν = 2 and 4.6 dB at ν ≈ 28.
\(\nu = 2K_{\mathrm{eff}}\)−5.7 / +16.0 dB at \(\nu = 2\)

\(\hat S/S \approx \chi^2_\nu/\nu\), so \(S \in [\nu\hat S/\chi^2_{\nu,0.975},\ \nu\hat S/\chi^2_{\nu,0.025}]\) (27.8). For \(\nu = 2\) the quantiles are 0.0506 and 7.378, giving −5.67 and +15.97 dB. The width depends only on \(\nu\), hence constant on a dB axis. Coverage measured 0.947 (\(\nu = 2\)) and 0.941 (\(\nu = 28.5\)).

Approx. modelχ² quantiles from the Wilson–Hilferty approximation (exact at ν = 2; error 0.3 dB at ν = 4, about 0.01 dB at ν = 16 and smaller beyond).
\(\nu\) is an equivalent count; correlated looks (overlapping lines, adjacent beams) make the nominal \(\nu\) optimistic.
A log of an average is biased low by ψ(K) − ln K; it cancels in a ratio only if both spectra had the same number of independent looks.
−2.51 dB at \(K = 1\)−0.07 dB at \(K = 32\)var = \(\psi_1(K) \approx 1/K\)

\(\ln\) of the mean of \(K\) unit exponentials has mean \(\psi(K) - \ln K\) and variance \(\psi_1(K)\). In an SLD both gates carry the same offset and it drops out. If the distal gate is partly noise, or its lines are more correlated (wider beam out of focus), the offsets differ and appear as a depth-dependent bias.

Exact formula
Exact for independent exponential bins; for heavily averaged spectra (\(K \gtrsim 30\)) the bias is below 0.1 dB and is dominated by other errors.
Multitaper: about 2NW − 1 orthogonal Slepian tapers, each concentrated in ±W; variance 1/K, resolution 2W.
\(N = 128,\ NW = 4\): \(K = 7\)\(2W\) = 2.5 MHz\(\lambda_6 = 0.937,\ \lambda_7 = 0.699\)

The DPSS are the eigenvectors of the sinc concentration kernel (27.9), computed stably from a commuting tridiagonal matrix (27.10). The first \(2NW - 1\) have \(\lambda \approx 1\) (little leakage); beyond that \(\lambda\) collapses. Averaging \(K\) eigenspectra gives \(\nu \approx 2K\) while using every sample of the gate.

SimulationEigenvalues of the 128 × 128 tridiagonal matrix (implicit QL), concentrations checked against the sinc kernel, JS.
For a smooth speckle spectrum, variance per unit bandwidth equals Hann–Welch's; the gain is leakage control and full use of short gates. \(2W\) smoothing biases peaked spectra.
Use Burg AR when the gate is shorter than the structure you need to resolve, and never read power from AR peak heights.
0.4 MHz in 64 samplesBurg P = 0.66 at \(p = 16\)Yule–Walker 0.09

Fourier resolution is \(f_s/N = 0.625\) MHz for \(N = 64\) at 40 MHz. Burg fits reflection coefficients to forward and backward errors without windowing (27.14), so it resolves below \(1/T\) at good SNR; Yule–Walker's biased autocorrelation is an implicit triangular window and smears. AIC picked \(p = 16\) (median).

SimulationOne simulated 64-sample gate (two tones, 20 dB SNR): Burg AR(16) and Hann periodogram, JS.
Resolution degrades fast with SNR; high orders give spurious peaks; Burg shows line splitting and phase-dependent bias at high SNR.
The log spectral ratio between depths is −4α(f)Δz: two for the round trip, two for power, α in nepers.
1 Np = 8.686 dB0.5 dB/cm/MHz = 0.0576 Np/cm/MHz−10 dB at 5 MHz, 2 cm

Amplitude decays as \(e^{-\alpha f z}\) one way, so power over the round trip decays as \(e^{-4\alpha f z}\) (27.21). 0.0576 × 4 × 5 × 2 = 2.30 Np of ln power = 10.0 dB. In dB directly: \(2\alpha_{\mathrm{dB}} f\Delta z\) = 2 × 0.5 × 5 × 2 = 10 dB. m27_attenuation_sld.m recovers the factor 4 by simulation (ratio 1.008).

Approx. modelLines: exact ln ratios. Dots: Gaussian scatter with the ln-spectrum std \(\sqrt{2\psi_1(32)}\) of a 32-look difference.
Assumes \(\alpha\) linear in \(f\) inside the band (\(\alpha \propto f^y\) with \(y\) ≈ 1–1.2 in tissue); with amplitude spectra the factor is 2.
Divide by a reference phantom: system response, focusing, gain and TGC cancel; attenuation and BSC differences remain.
same probe, focus, gain, TGCsame \(c\)

\(S(f,z) = G\cdot D\cdot e^{-4\int\alpha}\cdot\mathrm{BSC}\cdot C\) (27.16). Every factor that depends only on the system and geometry divides out (27.17) (Yao, Zagzebski & Madsen, 1990). In the toy focused simulation, the regression gives 0.70 dB/cm/MHz with the reference and 0.96 without it.

Approx. modelToy focal gain D = 1/(1 + ((z − z_F)/z_R)²), Gaussian 5 MHz 60 % spectrum, BSC ∝ f⁴; expected spectra, no noise.
Breaks when sample and reference sound speeds differ (different focus), when the reference \(\alpha\) or BSC is mis-specified, or when overlying layers differ (for BSC; the depth slope within a homogeneous ROI is immune to overlying attenuation).
Two-gate SLD error scales as 1/(Δz √(K T B³)).
\(\sigma \approx \sqrt{1.5\gamma\psi_1(K)/(B^3T)}/\Delta z\)0.069 dB/cm/MHz

Each ln spectrum has variance \(\psi_1(K)\) per independent frequency sample, \(BT/\gamma_w\) such samples in the band, and a least-squares slope over \(B\) has variance \(12\sigma^2/(nB^2)\). Divide by \(4\Delta z\) (27.22). \(K = 32\), \(T = 5.2\) µs, \(B = 4\) MHz, \(\Delta z = 2\) cm, Hann (\(\gamma = 1.75\)): 0.069 predicted, 0.067 measured.

Approx. modelVariance formula (27.22) with γ_w = 1.75 from Monte Carlo; independent looks, no electronic noise.
\(\gamma_w\) is a Monte Carlo constant for the window; additive noise, line correlation and a non-flat ln-spectrum variance across the band all raise \(\sigma\).
For a fixed ROI, precision is set by ROI length, band and lines: \(\sigma \propto (K B^3 L^3)^{-1/2}\), and the gate length drops out.
\(\sigma = \sqrt{4.5\gamma\psi_1(K)c/(B^3L^3)}\)\(L\) = 3 cm → 0.041 dB/cm/MHz

\(n = L/(cT/2)\) contiguous gates, each with frequency-slope variance \(\propto 1/T\), regressed over depth with \(\sum(z - \bar z)^2 \propto n^3(cT/2)^2\): \(T\) cancels (27.26). Doubling \(L\) gives 2.8×, doubling \(K\) 1.4×. Verified within 5 % by Monte Carlo for 2, 4 and 8 mm gates.

Approx. modelPrecision formula (27.26), Hann gates, K = 32, B = 4 MHz.
Needs a homogeneous ROI; overlapping tapered gates do somewhat better; very short gates violate the smooth-spectrum assumption and raise log-bias differences.
The centroid slope is −4α times the spectral variance at that depth, Gaussian or not.
\(\partial f_c/\partial z = -4\alpha_0\sigma^2(z)\)\(\sigma_p^2\) = 0.811 MHz² (5 MHz, 60 %)

\(\partial S/\partial z = -4\alpha_0 f S\), so \(\partial f_c/\partial z = -4\alpha_0(\langle f^2\rangle - \langle f\rangle^2)\) (27.18). Only a Gaussian keeps \(\sigma^2\) constant and the downshift linear. Using the transmit \(\sigma_p^2\) when the BSC (\(f^n\)), diffraction or the transducer reshape the spectrum biases \(\alpha\) in proportion to the variance error.

Exact formulaNumerical moments of the closed-form spectra.
Requires \(\alpha \propto f\); noise floors and band truncation of the centroid integral distort \(\sigma^2\) further.
Centroid and log-ratio methods see the same physics: the centroid trades the reference for a spectral-shape model.
both std 0.030 dB/cm/MHzcentroid bias +6 %

In m27_attenuation_sld.m, over the same 31 gates and 32 lines, the log-ratio regression gave 0.702 ± 0.030 and the reference-subtracted centroid 0.746 ± 0.030 dB/cm/MHz (true 0.70): equal precision, but the centroid kept a residual bias from focusing interacting with spectral shape. Raw centroid without a reference: 1.13.

SimulationMeans ± std reported by m27_attenuation_sld.m (Octave, toy focusing), drawn as Gaussians.
Numbers are for a toy focusing model; the centroid is attractive when no reference exists, if the spectral shape is known and stable.
BSC is in sr⁻¹cm⁻¹ and rises as f⁴ for small scatterers; the form factor bends it near ka ≈ 1.
\(k^4F\) peaks at \(ka = 1.555\)\(n = 4 - 1.654(ka)^2\)

Gaussian form factor \(F = \exp(-0.827k^2a^2)\): \(d\ln(k^4F)/d\ln f = 4 - 2\cdot 0.827(ka)^2\) (27.25). ESD 100 µm at 5 MHz: \(ka = 1.02\), \(n = 2.3\); ESD 200 µm: \(n = -2.9\) (the BSC already falls with frequency).

Exact formula
Other microstructures need other form factors (fluid sphere, exponential, spherical shell); the Gaussian is a convenient fit, not a law.
ESD is identifiable only when the band reaches ka of roughly 0.5 to 1; below that the fit returns noise.
ESD 20 µm, 3–7 MHz: \(ka \le 0.29\)36 % fits fail

For \(ka \ll 1\) the form factor is flat to within the spectral noise, so the \(k^2\) slope is poorly determined and its sign is random. m27_bsc_esd_fit.m (\(K = 32\), 4 mm gate): ESD 20 µm gives 36 % failures and a biased 33 ± 13 µm; ESD ≥ 60 µm (\(ka_{\max} \ge 0.85\)) is unbiased with std 11 µm down to 2 µm as \(ka\) grows.

Exact formulaGaussian form factor in dB; grey band = ln-spectrum std at K = 32 from the trigamma function.
Higher frequency extends \(ka\) but costs penetration and SNR; the usable band ends where the noise floor starts.
Attenuation and size trade off: an attenuation-compensation error tilts the spectrum and moves the ESD.
+0.1 dB/cm/MHz at 3 cmESD 100 → 78 µm

A compensation error \(\delta\alpha\) multiplies the BSC ratio by \(\exp(4\delta\alpha f z)\), a linear-in-\(f\) term that the \(k^2\) regression partly absorbs. With \(\delta\alpha\) = +0.1 dB/cm/MHz at \(z\) = 3 cm (≈ 0.14 Np per MHz of band) a true 100 µm reads 77.6 µm and 200 µm reads 190 µm (noise-free).

Exact formulaNoise-free least-squares fit of (27.24) with the tilt \(\exp(4\delta\alpha f z)\).
Smaller at shallow depth and for large \(ka\); joint estimation of \(\alpha\) and ESD from the same data reduces but correlates the errors.
Midband fit is the most precise Lizzi parameter because a least-squares line is best determined at the mean abscissa.
var at \(\bar f\): \(\sigma^2/n\)slope and midband uncorrelated

For a line fitted to \(n\) points, \(\mathrm{var}(\hat y(f)) = \sigma^2[1/n + (f - \bar f)^2/S_{ff}]\). At band center the second term vanishes and the estimate is uncorrelated with the slope; the intercept at 0 MHz is an extrapolation with large variance, which is why intercept images look noisier than midband images.

Exact formula
Lizzi parameters need the same attenuation compensation as the BSC; their interpretation in terms of size and concentration assumes a scattering model.

27.11Cheat sheet

ConceptFormulaUltrasound use
Windowed periodogram\(P = \lvert\sum w x e^{-j2\pi fn}\rvert^2/\sum w^2\); \(\mathrm{E}\{P\} = S * H_w\)Gate spectrum; leakage sets the usable dynamic range
Periodogram statistics\(P \approx S\chi^2_2/2\); std/mean = 1; dB: −2.51 ± 5.57Why one A-line's spectrum is useless alone
Welch variance\(1/K_{\mathrm{eff}} = (1/N_b)[1 + 2\sum(1 - m/N_b)\rho^2(mD)]\); Hann 50 %: \(\rho = 1/6\)Depth sub-gates or lateral lines as looks
\(\chi^2\) interval\([\nu\hat S/\chi^2_{\nu,.975},\ \nu\hat S/\chi^2_{\nu,.025}]\), \(\nu = 2K_{\mathrm{eff}}\)Error bars on spectra and BSC curves
Log of a K-look averagemean \(\psi(K) - \ln K\); var \(\psi_1(K) \approx 1/K\)Noise of every log-spectral QUS method
MultitaperDPSS of length \(N\), half-bandwidth \(W\); \(K \approx 2NW - 1\); var ≈ \(1/K\); resolution \(2W\)Short gates, high dynamic range spectra
AR / Burg\(S = \sigma^2/\lvert 1 + \sum a_k e^{-j2\pi fk}\rvert^2\); \(\mathrm{AIC} = N\ln\sigma_p^2 + 2p\)Features and moments from 1–2 mm gates
QUS signal model\(S = G\cdot D\cdot e^{-4\int\alpha\, dz}\cdot\mathrm{BSC}\cdot C\)What a reference must cancel
Log spectral ratio\(\ln(S_s/S_r) = \ln(\mathrm{BSC}_s/\mathrm{BSC}_r) - 4\Delta\alpha(f)z\) (Np, power)SLD and reference phantom attenuation
Units1 Np = 8.686 dB; ln P-ratio \(-4\alpha_{\mathrm{Np}}f\Delta z\) = \(-2\alpha_{\mathrm{dB}}f\Delta z\) dBConverting slopes to dB/cm/MHz
Two-gate SLD precision\(\sigma \approx \sqrt{1.5\gamma_w\psi_1(K)/(B^3T)}/\Delta z\); \(\gamma\) = 0.93 rect, 1.75 HannSizing \(K\), gate and \(\Delta z\)
ROI regression precision\(\sigma \approx \sqrt{4.5\gamma_w\psi_1(K)c/(B^3L^3)}\), independent of \(T\)Attenuation-map resolution vs precision
Centroid identity\(\partial f_c/\partial z = -4\alpha_0\sigma^2(z)\); Gaussian: \(f_c = f_0 - 4\alpha_0\sigma_p^2 z\)Spectral-shift attenuation
BSC from reference\(\mathrm{BSC}_s = \mathrm{BSC}_r(S_s/S_r)e^{4(\alpha_s - \alpha_r)z}\)Absolute backscatter, sr⁻¹cm⁻¹
Gaussian form factor\(F = \exp(-0.827k^2a^2)\); ESD = \(2a\); \(n = 4 - 1.654(ka)^2\)Effective scatterer size
Lizzi linedB spectrum ≈ slope·\(f\) + intercept; midband = value at band centerSpectral parameter imaging

27.12Interview and qualifying-exam questions

I1Why doesn't a longer gate make the periodogram less noisy, and what does? · screen

Each DFT bin of Gaussian data is circular complex Gaussian, so \(\lvert X\rvert^2\) is exponential with std equal to its mean. A longer record adds bins at spacing \(1/T\); it does not average within a bin. Variance falls only by averaging independent looks: segments (Bartlett/Welch), tapers (multitaper), adjacent bins (smoothing), or in QUS independent A-lines. Each costs resolution or spatial extent. Follow-up: "how many independent looks does a Welch estimate with Hann 50 % overlap have?" \(K_{\mathrm{eff}}\) from the overlap formula (27.7), e.g. 14.3 for 15 segments.

I2Derive the factor in the log spectral ratio and convert a measured slope to dB/cm/MHz. · screen

Amplitude decays as \(e^{-\alpha_0 f z}\) one way (\(\alpha_0\) in Np/cm/MHz); the echo travels \(2z\), and power squares it: \(e^{-4\alpha_0 f z}\). Between gates \(\Delta z\) apart, \(\ln S(z_2) - \ln S(z_1) = -4\alpha_0 f\Delta z\). If that difference has slope −0.46 per MHz for \(\Delta z\) = 2 cm: \(\alpha_0\) = 0.46/8 = 0.0575 Np/cm/MHz = 0.50 dB/cm/MHz. Follow-up: "and with amplitude spectra?" The slope is −0.23 per MHz and you divide by \(2\Delta z\): same answer.

I3What does a reference phantom cancel, and what does it not? · onsite

It cancels the transducer response, focusing/diffraction, gain, TGC and gate/window factors, because they multiply both spectra identically at the same depth (27.17). It does not cancel a sound-speed mismatch (the focus and beam shape move), differences in overlying attenuation (critical for BSC, irrelevant to the depth slope inside a homogeneous ROI), errors in the phantom's own \(\alpha_r\) and \(\mathrm{BSC}_r\) (they pass straight into the result), nonlinear propagation differences at high pressure, or electronic noise. Follow-up: "how would you check?" Scan a second phantom of known properties as a sample, and check the ln ratio for curvature in depth.

I4You need σ(α) ≤ 0.05 dB/cm/MHz with a 3–7 MHz band. Size the ROI. · onsite

Use \(\sigma \approx \sqrt{4.5\gamma\psi_1(K)c/(B^3L^3)}\) in Np (27.26), Hann \(\gamma = 1.75\), \(c\) = 0.154 cm/µs, \(B\) = 4 MHz. Target 0.05/8.686 = 0.00576 Np. With \(K = 32\) (\(\psi_1\) = 0.0317): \(L^3\) = 4.5·1.75·0.0317·0.154/(64·0.00576²) ≈ 18.2 → \(L\) ≈ 2.6 cm. With \(K = 16\), \(L\) ≈ 3.3 cm. Lines must be independent, i.e. spaced by about a beam width (\(\lambda F_\#\) ≈ 0.6 mm at 5 MHz, F/2): 32 independent lines span about 2 cm laterally; denser lines are correlated and count for less than their number. Follow-up: "what if the band shrinks to 2 MHz at depth?" \(\sigma\) grows by \(2^{1.5}\) = 2.8×, so \(L\) must grow by 2.

I5A 2 mm gate: Welch, multitaper or AR? · onsite

A 2 mm gate is about 104 samples at 40 MHz, so Fourier resolution is ≈ 0.4 MHz. For calibrated power (BSC, SLD) use multitaper (\(NW\) ≈ 3–4, \(K\) = 5–7): full use of the gate, \(\nu \approx 2K\), controlled leakage; Welch would split an already short gate. For locating features below \(1/T\) or for moments, Burg AR with AIC-chosen order; never use AR peak heights as power. In all cases the dominant averaging still comes from lateral lines. Follow-up: "what does \(NW\) cost?" Spectral smoothing over \(2W = 2NW f_s/N\).

I6Centroid downshift or spectral log difference? · onsite

\(\partial f_c/\partial z = -4\alpha_0\sigma^2(z)\) holds exactly (27.18), so the centroid needs the spectral variance at every depth; using the transmit \(\sigma_p^2\) assumes a Gaussian, unchanged by BSC, diffraction and noise. The log-ratio needs a reference but no shape model. On the same data both reached 0.030 dB/cm/MHz std; the centroid kept a 6 % bias. Choose the centroid when no reference exists and the shape is stable (or estimate \(\sigma^2\) per depth), the log ratio otherwise. Follow-up: "how does the BSC frequency dependence affect each?" It cancels in the SLD depth slope, but changes \(\sigma^2\) and biases the centroid with a fixed \(\sigma_p^2\).

I7Why does electronic noise bias attenuation estimates low, and how do you guard against it? · research

Noise adds a roughly flat power floor. Deep gates and high frequencies, where the signal is weakest, are lifted the most, so the distal spectrum looks less attenuated and the ln-ratio slope shrinks. Guards: measure the noise spectrum (no transmit) and subtract it in power before taking logs; restrict the band per depth to SNR above ~10 dB; weight the regression by the estimated variance; report the usable depth. Follow-up: "does noise subtraction bias anything?" It can produce negative powers at low SNR; clip or drop those bins, which reintroduces a small bias.

I8When can you estimate an effective scatterer diameter, and what correlates with it? · research

When the band covers \(ka\) of roughly 0.5 to 1.5, where the Gaussian form factor bends the \(f^4\) trend (local exponent \(4 - 1.654(ka)^2\), (27.25)). Below \(ka\) ≈ 0.3 the fit is noise (36 % failures for 20 µm at 3–7 MHz). The estimate correlates with the attenuation compensation (an error of 0.1 dB/cm/MHz at 3 cm takes 100 µm to 78 µm) and with the assumed form factor. Follow-up: "how would you pick the transducer?" Put the band so \(ka\) spans the size range of interest, then check SNR at the needed depth.

I9How would you turn SLD into an attenuation image with better resolution? · research

Write the log-spectral data from many small gates as a linear model \(y = A\alpha + b + \varepsilon\), with \(\alpha\) a spatial map (and optionally BSC-ratio terms as nuisance parameters), and solve it with regularization: Tikhonov or total variation in space, the same machinery as Module 13. Regularized SLD (Coila & Lavarello, IEEE TUFFC, 2018) is one published example. Report bias and resolution of the regularized map (model-resolution matrix), not just its smoothness. Follow-up: "what is the price?" Bias toward the prior: TV shrinks contrast of small inclusions, Tikhonov blurs edges.

I10Why might a log-averaged spectrum differ systematically between two depths even with no attenuation difference? · onsite

The expected log of an average of \(K\) exponentials, \(\psi(K) - \ln K\), depends on the effective number of independent looks. If the distal lines are more correlated (beam widens past focus) or contain noise, \(K_{\mathrm{eff}}\) differs between depths and the offsets no longer cancel: −2.51 dB at \(K = 1\) vs −0.07 dB at \(K = 32\). Fix by averaging power before the log, matching \(K_{\mathrm{eff}}\), or correcting with \(\psi(K_{\mathrm{eff}})\). Follow-up: "average log or log average?" Log of the average: lower variance and a bias that vanishes as \(1/(2K)\).

I11How would you validate a new attenuation estimator before a 510(k) or a paper? · research

Phantoms with \(\alpha\) measured independently by insertion loss (narrowband substitution) over the same band; several \(\alpha\) values and depths; bias and standard deviation versus ROI size, not only means; a layered phantom to test overlying-layer immunity; a sound-speed mismatch case; and in vivo repeatability (same subject, repeated scans, different operators). Report units, band, ROI dimensions and the number of independent lines. Public benchmarks that pair attenuation with sound-speed ground truth help comparison across methods. Follow-up: "what accuracy is enough?" Set it from the clinical decision threshold and the population spread, then size the ROI with the precision formula.

27.13MATLAB 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: m27_periodogram_welch.m · m27_multitaper_dpss.m · m27_ar_burg.m · m27_attenuation_sld.m · m27_bsc_esd_fit.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m27_periodogram_welch.mPeriodogram χ²₂ statistics, Welch variance formula, χ² confidence intervals and their coverage · change N, M · runtime ≈ 3 s 
%% m27_periodogram_welch.m  -  periodogram bias/variance, chi-square statistics, Welch averaging
% Purpose : speckle-like RF (white Gaussian scatterers convolved with a 5 MHz,
%           60 % bandwidth Gaussian pulse, fs = 40 MHz) has a known power
%           spectral density S(f). Show that
%           (1) the periodogram P(f) is unbiased for smooth S but its relative
%               standard deviation stays 1 however long the record (not
%               consistent); P/S is exponential (chi-square with 2 dof / 2);
%           (2) on a dB scale one bin has std 5.57 dB and sits 2.51 dB low;
%           (3) Welch averaging (Hann, 50 % overlap) cuts the variance to
%               (1 + 2 sum (1-m/Nb) rho(m)^2)/Nb, rho = taper overlap correlation,
%               while the bin width grows from fs/N to fs/M;
%           (4) chi-square confidence intervals with nu = 2 K_eff cover the truth
%               95 % of the time.
% Change  : N (record length), M (Welch segment length), number of trials.
% Expect  : std(P/S) about 1.00 for N = 64, 256, 1024; dB mean -2.5, std 5.6;
%           Welch N = 1024, M = 128: Nb = 15, variance 0.070 predicted, 0.074 measured;
%           CI coverage about 0.95 for both.
clear; close all; rng(1);

%% Parameters
fs = 40;  f0 = 5;  bw = 0.6;                 % MHz; -6 dB fractional bandwidth (amplitude)
sigf = bw*f0/(2*sqrt(2*log(2)));             % std of the amplitude spectrum (MHz)
Nfft = 4096;                                 % length of the generated stationary record
fF = (0:Nfft-1)'/Nfft*fs;  fF(fF >= fs/2) = fF(fF >= fs/2) - fs;
H = exp(-(abs(fF) - f0).^2/(2*sigf^2));      % amplitude transfer function (two-sided, real, even)
Strue = @(f) exp(-(abs(f) - f0).^2/sigf^2);  % PSD of x for unit-variance white input (per sample, fs = 1 units)
% x = ifft(fft(white).*H) has PSD |H|^2 in "per bin" units; P = |X|^2/N estimates it.
makeX = @(n) real(ifft(fft(randn(Nfft, n)).*repmat(H, 1, n)));
inband = @(f) abs(f - f0) <= 1.5;            % +/-1.5 MHz around f0 (inside the -14 dB band)

%% 1. Periodogram: variance does not fall with N
ntr = 200;
fprintf('Periodogram (rectangular window), %d realizations, bins within f0 +/- 1.5 MHz:\n', ntr);
fprintf('    N    bin width (MHz)   mean(P/S)   std(P/S)   mean dB   std dB\n');
for N = [64 256 1024]
  X = makeX(ntr);  x = X(1:N, :);
  P = abs(fft(x)).^2/N;                      % periodogram, per-bin units
  f = (0:N-1)'/N*fs;  k = inband(f) & f < fs/2;
  R = P(k, :)./repmat(Strue(f(k)), 1, ntr);  % should be Exp(1)
  RdB = 10*log10(R(:));
  fprintf('%5d   %10.3f       %8.3f   %8.3f   %7.2f   %6.2f\n', N, fs/N, mean(R(:)), std(R(:)), mean(RdB), std(RdB));
  if N == 1024, Rkeep = R(:); end
end
fprintf('Theory: mean 1, std 1; dB mean = -gamma*10/ln10 = %.2f dB, dB std = (pi/sqrt 6)*10/ln10 = %.2f dB\n', ...
        -0.5772156649*10/log(10), pi/sqrt(6)*10/log(10));

%% 2. Welch: Hann segments of length M, 50 % overlap
N = 1024;  M = 128;  S = M/2;  starts = 1:S:(N - M + 1);  Nb = numel(starts);
w = 0.5 - 0.5*cos(2*pi*((0:M-1)' + 0.5)/M);  U = sum(w.^2);
rho = zeros(Nb-1, 1);
for m = 1:Nb-1
  sh = m*S;  if sh < M, rho(m) = sum(w(1:M-sh).*w(1+sh:M))/U; end
end
varTheory = (1 + 2*sum((1 - (1:Nb-1)'/Nb).*rho.^2))/Nb;
X = makeX(ntr);  fW = (0:M-1)'/M*fs;  kW = inband(fW) & fW < fs/2;
Pw = zeros(M, ntr);
for b = 1:Nb
  seg = X(starts(b):starts(b)+M-1, :).*repmat(w, 1, ntr);
  Pw = Pw + abs(fft(seg)).^2/U;
end
Pw = Pw/Nb;
% expected value = S convolved with the Hann spectral window; S is smooth, so compare with S directly
Rw = Pw(kW, :)./repmat(Strue(fW(kW)), 1, ntr);
fprintf('\nWelch N = %d, M = %d (Hann, 50%% overlap): Nb = %d segments, rho(M/2) = %.4f\n', N, M, Nb, rho(1));
fprintf('  bin width fs/M = %.3f MHz (vs %.3f MHz raw)\n', fs/M, fs/N);
fprintf('  variance of P/S: measured %.4f, theory %.4f, K_eff = %.1f, 1/Nb = %.4f\n', var(Rw(:)), varTheory, 1/varTheory, 1/Nb);
fprintf('  mean of P/S = %.3f (window-smoothing bias is small for this smooth spectrum)\n', mean(Rw(:)));

%% 3. Chi-square confidence intervals and their coverage
% 95 % quantiles of chi2 with nu dof via the Wilson-Hilferty approximation (base MATLAB has no chi2inv)
z975 = 1.959964;
chi2q = @(nu, z) nu*(1 - 2/(9*nu) + z*sqrt(2/(9*nu))).^3;
% exact for nu = 2: chi2_2 quantile = -2 ln(1-p)
lo2 = -2*log(1 - 0.025);  hi2 = -2*log(1 - 0.975);
cov2 = mean(2*Rkeep >= lo2 & 2*Rkeep <= hi2);
fprintf('\nRaw periodogram (nu = 2): 95%% interval for S is [P*2/%.3f, P*2/%.4f] = [%.2f dB, +%.2f dB] around P\n', ...
        hi2, lo2, 10*log10(2/hi2), 10*log10(2/lo2));
fprintf('  span %.1f dB; empirical coverage %.3f\n', 10*log10(hi2/lo2), cov2);
nu = 2/varTheory;  lo = chi2q(nu, -z975);  hi = chi2q(nu, z975);
covW = mean(nu*Rw(:) >= lo & nu*Rw(:) <= hi);
fprintf('Welch (nu = 2 K_eff = %.1f): interval [%.2f dB, +%.2f dB], span %.1f dB; coverage %.3f\n', ...
        nu, 10*log10(nu/hi), 10*log10(nu/lo), 10*log10(hi/lo), covW);

%% Figures
x = X(:, 1);  P1 = abs(fft(x(1:N))).^2/N;  f1 = (0:N-1)'/N*fs;
Pw1 = zeros(M, 1);
for b = 1:Nb, Pw1 = Pw1 + abs(fft(x(starts(b):starts(b)+M-1).*w)).^2/U; end
Pw1 = Pw1/Nb;
figure('Position', [100 100 900 360]);
subplot(1, 2, 1);
h1 = f1 < fs/2;  h2 = fW < fs/2;  fT = linspace(0, fs/2, 800)';
plot(f1(h1), 10*log10(P1(h1)), 'Color', [0.6 0.6 0.6]); hold on;
plot(fW(h2), 10*log10(Pw1(h2)), 'b-o', 'MarkerSize', 3, 'LineWidth', 1.2);
plot(fT, 10*log10(Strue(fT)), 'r--', 'LineWidth', 1.5);
xlim([0 12]); ylim([-40 12]); grid on;
xlabel('frequency (MHz)'); ylabel('PSD (dB re peak)');
legend('periodogram N = 1024', 'Welch M = 128, 50 %', 'true S(f)', 'Location', 'northeast');
title('One record: raw periodogram vs Welch');
subplot(1, 2, 2);
edges = 0:0.25:6;  cnt = histc(Rkeep, edges);  cnt = cnt(1:end-1);
bar(edges(1:end-1) + 0.125, cnt/(numel(Rkeep)*0.25), 1, 'FaceColor', [0.7 0.8 1]); hold on;
xx = linspace(0, 6, 200);  plot(xx, exp(-xx), 'r', 'LineWidth', 1.5);
xlabel('P(f) / S(f)'); ylabel('probability density'); grid on;
title('Periodogram bins are exponential (\chi^2_2/2)');
legend('N = 1024, in-band bins', 'e^{-x}');

What to look for: \(\mathrm{std}(P/S)\) ≈ 1.00 for \(N\) = 64, 256, 1024; dB mean −2.6 (theory −2.51) and std 5.6 (5.57); Welch \(N_b\) = 15, variance 0.074 vs 0.070 predicted; interval spans 21.6 dB (\(\nu = 2\)) and 4.6 dB (\(\nu = 28.5\)) with coverage 0.947 and 0.941.

m27_multitaper_dpss.mDPSS from the tridiagonal matrix via eig, concentrations against the sinc kernel, multitaper vs Hann vs Welch on a speckle gate · change N, NW · runtime ≈ 3 s 
%% m27_multitaper_dpss.m  -  Slepian (DPSS) tapers from eig, and Thomson's multitaper spectrum
% Purpose : (1) build the K = 2NW - 1 discrete prolate spheroidal sequences of
%           length N without a toolbox, as eigenvectors of the symmetric
%           tridiagonal matrix that commutes with the sinc (prolate) kernel;
%           (2) check them against the sinc-kernel eigenproblem: the
%           concentration lambda_k (fraction of taper energy inside |f| <= W)
%           equals the sinc-kernel eigenvalue;
%           (3) estimate the spectrum of a short speckle gate with K tapers and
%           verify that the variance falls as 1/K (nu = 2K dof);
%           (4) compare variance x bandwidth with Hann-Welch at 50 % overlap.
% Change  : N, NW (half bandwidth in bins), K.
% Expect  : NW = 4, N = 128: lambda = 1.0000 ... 0.9994 0.9926 0.9369, then
%           lambda_7 = 0.70 (leaky); orthonormality error ~1e-15;
%           var(S_mt/S) = 0.17, a little above 1/K = 0.143 because the speckle
%           spectrum is not flat across 2W = 2NW*fs/N = 2.5 MHz; the same
%           smoothing biases the mean low (0.91) on this short gate. Hann-Welch
%           with the same spectral-window width gives about the same variance.
clear; close all; rng(1);

%% Parameters
N  = 128;            % gate length (samples); 3.2 us at 40 MHz = 2.46 mm of depth
NW = 4;              % time-half-bandwidth product
W  = NW/N;           % half bandwidth (cycles/sample)
K  = 2*NW - 1;       % number of well-concentrated tapers
fs = 40;             % MHz

%% 1. DPSS from the tridiagonal matrix (Slepian 1978; Percival & Walden 1993, ch. 8)
n = (0:N-1)';
d = ((N - 1)/2 - n).^2*cos(2*pi*W);          % diagonal
e = n(2:end).*(N - n(2:end))/2;              % off-diagonal
T = diag(d) + diag(e, 1) + diag(e, -1);
[V, D] = eig(T);  [~, idx] = sort(diag(D), 'descend');
V = V(:, idx(1:K+1));                        % keep one extra taper to show the leakage cliff
for k = 1:K+1                                % sign convention: even tapers sum > 0, odd tapers start > 0 slope
  if mod(k-1, 2) == 0, s = sign(sum(V(:, k))); else, s = sign(sum((n - (N-1)/2).*V(:, k))); end
  V(:, k) = s*V(:, k);
end
fprintf('Orthonormality: max |V''V - I| = %.1e\n', max(max(abs(V'*V - eye(K+1)))));

%% 2. Concentration from the sinc kernel  A(m,n) = sin(2 pi W (m-n)) / (pi (m-n))
[mm, nn] = meshgrid(n, n);  dmn = mm - nn;
A = sin(2*pi*W*dmn)./(pi*dmn + (dmn == 0));  A(dmn == 0) = 2*W;
lam = diag(V'*A*V);                          % Rayleigh quotient = energy fraction in |f| <= W
lamA = sort(eig((A + A')/2), 'descend');
fprintf('\n k   lambda_k (from tapers)   eig of sinc kernel   out-of-band energy (dB)\n');
for k = 1:K+1
  fprintf('%2d     %.6f             %.6f            %7.1f\n', k-1, lam(k), lamA(k), 10*log10(max(1 - lam(k), 1e-300)));
end
fprintf('Number of eigenvalues > 0.5: %d (Shannon number 2NW = %d)\n', sum(lamA > 0.5), 2*NW);

%% 3. Multitaper spectrum of a speckle gate, Monte Carlo variance
f0 = 5;  bw = 0.6;  sigf = bw*f0/(2*sqrt(2*log(2)));
Nf = 1024;  ff = (0:Nf-1)'/Nf*fs;  ff(ff >= fs/2) = ff(ff >= fs/2) - fs;
H = exp(-(abs(ff) - f0).^2/(2*sigf^2));
Strue = @(f) exp(-(abs(f) - f0).^2/sigf^2);
ntr = 400;  Np = 512;  fP = (0:Np-1)'/Np*fs;  kb = abs(fP - f0) <= 1.2;
Rmt = zeros(sum(kb), ntr);  Rhw = Rmt;  Rp = Rmt;
w = 0.5 - 0.5*cos(2*pi*(n + 0.5)/N);  w = w/sqrt(sum(w.^2));
M = 32;  S = 16;  st = 1:S:(N - M + 1);  wh = 0.5 - 0.5*cos(2*pi*((0:M-1)' + 0.5)/M);  wh = wh/sqrt(sum(wh.^2));
for t = 1:ntr
  x = real(ifft(fft(randn(Nf, 1)).*H));  x = x(1:N);
  Smt = mean(abs(fft(V(:, 1:K).*repmat(x, 1, K), Np)).^2, 2);
  Sp  = abs(fft(x.*w, Np)).^2;
  Sw  = zeros(Np, 1);
  for b = 1:numel(st), Sw = Sw + abs(fft(x(st(b):st(b)+M-1).*wh, Np)).^2; end
  Sw = Sw/numel(st);
  Rmt(:, t) = Smt(kb)./Strue(fP(kb));  Rp(:, t) = Sp(kb)./Strue(fP(kb));  Rhw(:, t) = Sw(kb)./Strue(fP(kb));
end
fprintf('\nSpeckle gate, N = %d (%.2f mm), bins within f0 +/- 1.2 MHz, %d trials:\n', N, N/fs*1.54/2, ntr);
fprintf('  Hann periodogram : var(S^/S) = %.3f, mean = %.3f\n', var(Rp(:)), mean(Rp(:)));
fprintf('  Multitaper K = %d: var(S^/S) = %.3f (1/K = %.3f), mean = %.3f\n', K, var(Rmt(:)), 1/K, mean(Rmt(:)));
fprintf('  Hann-Welch M = %d, 50%% overlap (%d segments): var = %.3f, mean = %.3f\n', M, numel(st), var(Rhw(:)), mean(Rhw(:)));
fprintf('  Resolution: multitaper 2W = %.3f MHz; Welch bin fs/M = %.3f MHz\n', 2*W*fs, fs/M);

%% 4. Bandwidth (int Wk)^2/int Wk^2 of the averaged spectral windows (in MHz)
L = 32*N;  fL = (0:L-1)'/L*fs;
Wmt = mean(abs(fft(V(:, 1:K), L)).^2, 2);  Wmt = Wmt/sum(Wmt);
Whw = abs(fft(wh, L)).^2;  Whw = Whw/sum(Whw);
Bmt = 1/sum(Wmt.^2)*fs/L;  Bhw = 1/sum(Whw.^2)*fs/L;
fprintf('  Spectral-window bandwidth: multitaper %.2f MHz, Welch segment %.2f MHz\n', Bmt, Bhw);
fprintf('  variance x bandwidth x T: multitaper %.2f, Welch %.2f (T = %.2f us)\n', ...
        var(Rmt(:))*Bmt*N/fs, var(Rhw(:))*Bhw*N/fs, N/fs);

%% Figures
figure('Position', [100 100 950 380]);
subplot(1, 2, 1);
plot(n, V(:, 1:4), 'LineWidth', 1.3); grid on; xlim([0 N-1]);
xlabel('sample n'); ylabel('taper value');
legend('v_0', 'v_1', 'v_2', 'v_3', 'Location', 'south');
title(sprintf('DPSS tapers, N = %d, NW = %g', N, NW));
subplot(1, 2, 2);
fW = fL;  fW(fW >= fs/2) = fW(fW >= fs/2) - fs;  [fW, o] = sort(fW);
Wk = abs(fft(V, L)).^2;  Wk = Wk(o, :);
plot(fW, 10*log10(Wk(:, 1)/max(Wk(:, 1))), 'b', 'LineWidth', 1.3); hold on;
plot(fW, 10*log10(Wk(:, K)/max(Wk(:, K))), 'm', 'LineWidth', 1.1);
plot(fW, 10*log10(Wk(:, K+1)/max(Wk(:, K+1))), 'Color', [0.6 0.6 0.6]);
plot([W W]*fs, [-120 5], 'k--', [-W -W]*fs, [-120 5], 'k--');
xlim([-6 6]); ylim([-120 5]); grid on;
xlabel('frequency offset (MHz)'); ylabel('|V_k(f)|^2 (dB)');
legend('k = 0', sprintf('k = %d', K-1), sprintf('k = %d (leaky)', K), '\pm W', 'Location', 'south');
title('Spectral windows: concentration in |f| \leq W');

What to look for: orthonormality error 1e−15; \(\lambda_{0..7}\) from both eigenproblems agree (… 0.9926, 0.9369, then 0.699); 8 eigenvalues above 0.5 (\(2NW = 8\)); \(\mathrm{var}(\hat S/S)\) = 1.08 (Hann periodogram), 0.17 (multitaper, \(1/K\) = 0.143), 0.17 (Welch); variance × bandwidth × \(T\) ≈ 1.4 for both.

m27_ar_burg.mBurg vs Yule–Walker vs periodogram resolution on a 64-sample gate, AIC order, AR centroid of speckle · change N, separation, SNR, p · runtime ≈ 8 s 
%% m27_ar_burg.m  -  parametric (AR) spectra for short gates: Burg vs Yule-Walker vs periodogram
% Purpose : (1) two tones 0.4 MHz apart in a 64-sample gate (1.6 us at 40 MHz,
%               1.23 mm of depth): the periodogram cannot separate them
%               (Rayleigh resolution fs/N = 0.625 MHz); Burg AR(p) can once p
%               is large enough; Yule-Walker (biased autocorrelation = implicit
%               triangular window) needs a much higher order;
%           (2) AIC order selection: N ln(E_p) + 2p;
%           (3) on a short speckle gate, the spectral centroid from Burg AR vs a
%               Hann periodogram: both unbiased, AR about 20 % lower std, and
%               the AR order hardly matters because the one-sided centroid is
%               fs/4 - (2 fs/pi^2) sum_{odd l} rho(l)/l^2, dominated by lag 1,
%               and every Burg model of order >= l shares rho(0..l).
% Change  : N, tone separation, SNR, model order p.
% Expect  : P(resolved) Burg 0.37 / 0.66 / 0.75 at p = 12 / 16 / 24;
%           Yule-Walker 0.00 / 0.09 / 0.29; Hann periodogram 0.02; AIC picks
%           p = 16 (median). Speckle centroid std 0.36 MHz (periodogram) vs
%           0.29 MHz (Burg, any order).
clear; close all; rng(1);  clear x1;

%% Parameters
fs = 40;  N = 64;  f1 = 4.8;  f2 = 5.2;     % MHz, samples
snr_dB = 20;                                % per tone: (A^2/2)/sigma^2 with A = 1
ntr = 200;  orders = [4 8 12 16 24];
fg = linspace(3, 7, 801)';                  % evaluation grid (MHz)
E1 = exp(-1i*2*pi*fg/fs*(0:max(orders)));   % e^{-j 2 pi f k / fs}
arpsd = @(a, E) E./abs(E1(:, 1:numel(a))*a).^2;
% Burg recursion is written out inline below (scripts cannot hold local functions portably)

%% 1. Resolution test
fprintf('Two tones %.1f and %.1f MHz, N = %d (%.2f us), SNR %d dB per tone, %d trials\n', f1, f2, N, N/fs, snr_dB, ntr);
fprintf('Periodogram Rayleigh resolution fs/N = %.3f MHz > separation %.1f MHz\n', fs/N, f2 - f1);
resB = zeros(numel(orders), 2);  resY = resB;  resP = 0;  aicSel = zeros(ntr, 1);
n = (0:N-1)';  sig = sqrt(0.5*10^(-snr_dB/10));
for t = 1:ntr
  x = cos(2*pi*f1/fs*n + 2*pi*rand) + cos(2*pi*f2/fs*n + 2*pi*rand) + sig*randn(N, 1);
  % --- Burg up to max order, storing every intermediate model
  pmax = max(orders);  ef = x;  eb = x;  a = 1;  Ep = zeros(pmax+1, 1);  Ep(1) = x'*x/N;  Aall = cell(pmax, 1);
  for m = 0:pmax-1
    ff = ef(m+2:N);  bb = eb(m+1:N-1);
    k = -2*(bb'*ff)/(ff'*ff + bb'*bb);
    a = [a; 0] + k*flipud([a; 0]);
    ef(m+2:N) = ff + k*bb;  eb(m+2:N) = bb + k*ff;
    Ep(m+2) = Ep(m+1)*(1 - k^2);  Aall{m+1} = a;
  end
  aic = N*log(Ep) + 2*(0:pmax)';  [~, im] = min(aic);  aicSel(t) = im - 1;
  % --- Yule-Walker with the biased autocorrelation estimate
  r = zeros(pmax+1, 1);  for l = 0:pmax, r(l+1) = x(1:N-l)'*x(1+l:N)/N; end
  for io = 1:numel(orders)
    p = orders(io);
    for meth = 1:2
      if meth == 1, ap = Aall{p}; Pp = arpsd(ap, Ep(p+1));
      else, ay = [1; -toeplitz(r(1:p))\r(2:p+1)]; Pp = arpsd(ay, r(1) + r(2:p+1)'*ay(2:end)); end
      pk = find(Pp(2:end-1) > Pp(1:end-2) & Pp(2:end-1) >= Pp(3:end) & Pp(2:end-1) > 0.01*max(Pp)) + 1;
      ok = numel(pk) == 2 && all(abs(sort(fg(pk)) - [f1; f2]) < 0.15);
      if meth == 1 && p == 16 && ok && ~exist('x1', 'var'), x1 = x; A1 = Aall; E1p = Ep; end
      if meth == 1, resB(io, :) = resB(io, :) + [ok, numel(pk) > 2]; else, resY(io, :) = resY(io, :) + [ok, numel(pk) > 2]; end
    end
  end
  Px = abs(exp(-1i*2*pi*fg/fs*n')*(x.*(0.5 - 0.5*cos(2*pi*(n + 0.5)/N)))).^2;
  pk = find(Px(2:end-1) > Px(1:end-2) & Px(2:end-1) >= Px(3:end) & Px(2:end-1) > 0.01*max(Px)) + 1;
  resP = resP + (numel(pk) == 2 && all(abs(sort(fg(pk)) - [f1; f2]) < 0.15));
end
fprintf('\n order p   Burg: P(resolved)  P(>2 peaks)   Yule-Walker: P(resolved)\n');
for io = 1:numel(orders)
  fprintf('   %3d        %5.3f          %5.3f              %5.3f\n', orders(io), resB(io, 1)/ntr, resB(io, 2)/ntr, resY(io, 1)/ntr);
end
fprintf(' Hann periodogram: P(resolved) = %.3f\n', resP/ntr);
qs = sort(aicSel);
fprintf(' AIC-selected Burg order: median %d, 10-90 %% range %d-%d\n', median(aicSel), qs(round(0.1*ntr)), qs(round(0.9*ntr)));

%% 3. Short speckle gate: centroid from Burg AR(p) vs Hann periodogram
f0 = 5;  sigf = 0.6*f0/(2*sqrt(2*log(2)));  Nf = 1024;
fF = (0:Nf-1)'*fs/Nf;  fF(fF >= fs/2) = fF(fF >= fs/2) - fs;
H = exp(-(abs(fF) - f0).^2/(2*sigf^2));
fc = linspace(0.5, 12, 1151)';  Ec = exp(-1i*2*pi*fc/fs*(0:16));
cen = @(P) sum(fc.*P)/sum(P);
ntr2 = 500;  Cp = zeros(ntr2, 1);  Cb = zeros(ntr2, 3);  pl = [4 8 16];
wN = 0.5 - 0.5*cos(2*pi*(n + 0.5)/N);  En = exp(-1i*2*pi*fc/fs*n');
for t = 1:ntr2
  x = real(ifft(fft(randn(Nf, 1)).*H));  x = x(1:N);
  Cp(t) = cen(abs(En*(x.*wN)).^2);
  ef = x;  eb = x;  a = 1;  Ep = x'*x/N;
  for m = 0:max(pl)-1
    ff = ef(m+2:N);  bb = eb(m+1:N-1);  k = -2*(bb'*ff)/(ff'*ff + bb'*bb);
    a = [a; 0] + k*flipud([a; 0]);  ef(m+2:N) = ff + k*bb;  eb(m+2:N) = bb + k*ff;  Ep = Ep*(1 - k^2);
    j = find(pl == m + 1);  if ~isempty(j), Cb(t, j) = cen(Ep./abs(Ec(:, 1:numel(a))*a).^2); end
  end
end
fprintf('\nSpeckle gate N = %d (%.2f mm), true power-spectrum centroid %.2f MHz, %d gates:\n', N, N/fs*1.54/2, f0, ntr2);
fprintf('  Hann periodogram : mean %.3f MHz, std %.3f MHz\n', mean(Cp), std(Cp));
for j = 1:numel(pl), fprintf('  Burg AR(%2d)      : mean %.4f MHz, std %.4f MHz\n', pl(j), mean(Cb(:, j)), std(Cb(:, j))); end

%% Figures
figure('Position', [100 100 950 380]);
subplot(1, 2, 1);
Px1 = abs(exp(-1i*2*pi*fg/fs*n')*(x1.*(0.5 - 0.5*cos(2*pi*(n + 0.5)/N)))).^2;
plot(fg, 10*log10(Px1/max(Px1)), 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2); hold on;
for p = [8 16]
  Pp = arpsd(A1{p}, E1p(p+1));  plot(fg, 10*log10(Pp/max(Pp)), 'LineWidth', 1.3);
end
plot([f1 f1], [-40 3], 'k:', [f2 f2], [-40 3], 'k:');
ylim([-40 3]); grid on; xlabel('frequency (MHz)'); ylabel('dB re peak');
legend('Hann periodogram', 'Burg p = 8', 'Burg p = 16', 'true tones', 'Location', 'south');
title(sprintf('A gate Burg(16) resolves: N = %d, tones %.1f MHz apart', N, f2 - f1));
subplot(1, 2, 2);
plot(orders, resB(:, 1)/ntr, 'o-', orders, resY(:, 1)/ntr, 's-', 'LineWidth', 1.3); grid on;
xlabel('AR order p'); ylabel('P(two peaks within 0.15 MHz)'); ylim([0 1]);
legend('Burg', 'Yule-Walker (biased r)', 'Location', 'northwest');
title('Resolution probability vs order');

What to look for: P(resolved) Burg 0.37 / 0.66 / 0.75 at \(p\) = 12 / 16 / 24, Yule–Walker 0.00 / 0.09 / 0.29, Hann periodogram 0.02; AIC median \(p = 16\); speckle centroid std 0.36 MHz (periodogram) vs 0.29 MHz (Burg, orders 4–16).

m27_attenuation_sld.mSpeckle A-lines with per-scatterer two-way attenuation and toy focusing; factor 4, reference method, two-gate SLD, centroid, and their spread · change αs, K, gate, band, focus · runtime ≈ 14 s 
%% m27_attenuation_sld.m  -  attenuation by the spectral log difference, with a reference phantom
% Purpose : simulate speckle RF A-lines in which EVERY scatterer's echo carries
%           its own two-way attenuation exp(-2 alpha f z) and a toy focusing
%           (diffraction) factor D(f,z), for a sample (alpha_s) and a reference
%           phantom (alpha_r) scanned with the same system. Then
%           (1) show that ln S_s(f,z) - ln S_r(f,z) falls with depth at slope
%               -4 (alpha_s - alpha_r) f  (POWER spectra, two-way, Np/cm);
%           (2) estimate alpha by the reference phantom method (regression over
%               depth, then over frequency) and by the two-gate spectral log
%               difference (SLD), and show the diffraction bias you get without
%               the reference;
%           (3) measure the estimator spread over independent groups of A-lines
%               and compare with sigma_alpha ~ sqrt(1.5 g_w psi1(K)/(B^3 T))/dz
%               (Np/cm/MHz; g_w = 1.75 for a Hann gate, 0.93 rectangular,
%               Monte Carlo on white noise; psi1 = trigamma, about 1/K);
%           (4) run the centroid-downshift estimator on the same data, raw
%               and with the reference phantom's centroid slope subtracted.
% Change  : alpha_s, gate length Lg, lines per estimate K, band, focus zF, snr_dB.
% Expect  : ratio of measured to -4 x dalpha slope = 1.00; reference method
%           0.70 dB/cm/MHz; without the reference the focusing factor biases
%           the estimate to about 0.96; two-gate SLD std 0.067 dB/cm/MHz for
%           K = 32, Lg = 4 mm, dz = 2 cm, B = 4 MHz (formula 0.069); regression
%           over 31 gates (1.5-4.5 cm) std 0.030; raw centroid 1.13 (focusing
%           bias), reference-subtracted centroid 0.75 (std 0.030, 6 % model bias).
clear; close all; rng(1);

%% Parameters
c   = 0.154;            % cm/us (1540 m/s)
fs  = 40;               % MHz
f0  = 5;  bw = 0.6;     % MHz, -6 dB fractional bandwidth of the AMPLITUDE spectrum
sigf = bw*f0/(2*sqrt(2*log(2)));   % amplitude-spectrum std (MHz); power std = sigf/sqrt(2)
NP  = 20/log(10);       % 8.686 dB per neper
alpha_s = 0.70;         % sample attenuation, dB/cm/MHz
alpha_r = 0.50;         % reference phantom attenuation, dB/cm/MHz (known)
zF  = 2.0;  Fnum = 2;   % focal depth (cm) and F-number of the toy focusing model
snr_dB = Inf;           % electronic SNR at the shallowest gate (Inf = no noise)
zmax = 5.0;             % deepest scatterer (cm)
Nlines = 4096;          % independent A-lines per medium
K   = 32;               % A-lines averaged per estimate (one ROI)
Lg  = 0.4;              % gate length (cm) = 4 mm
band = [3 7];           % analysis band (MHz)

%% Simulation: RF(f) = G(f) sum_n s_n sqrt(D(f,z_n)) exp(-2 alpha f z_n) exp(-j 2 pi f t_n)
Nt = 4096;  tn = (0:ceil(2*zmax/c*fs))'/fs;  zn = c*tn/2;      % scatterer times (us) and depths (cm)
k1 = round(0.5/fs*Nt);  k2 = round(10.5/fs*Nt);  fk = (k1:k2)'*fs/Nt;   % 0.5-10.5 MHz (G < -68 dB outside)
G  = exp(-(fk - f0).^2/(2*sigf^2));                            % amplitude transfer function
zR = 3.5*(c./fk)*Fnum^2;                                       % toy half depth of field 3.5 lambda F#^2 (cm), larger at low f
Dsq = 1./sqrt(1 + ((repmat(zn', numel(fk), 1) - zF)./repmat(zR, 1, numel(zn))).^2);  % sqrt of intensity factor
Ph  = exp(-1i*2*pi*fk*tn');                                    % propagation phase
tic;
Ms = repmat(G, 1, numel(zn)).*Dsq.*exp(-2*(alpha_s/NP)*fk*zn').*Ph;   % frequency x scatterer transfer, sample
Mr = repmat(G, 1, numel(zn)).*Dsq.*exp(-2*(alpha_r/NP)*fk*zn').*Ph;   % same for the reference phantom
makeRF = @(M, S) real(ifft([zeros(k1, size(S, 2)); M*S; ...
           zeros(Nt - k2 - 1, size(S, 2))]))*2;                % one-sided spectrum -> real RF
RFs = zeros(Nt, Nlines);  RFr = RFs;
for b = 1:Nlines/256
  idx = (b-1)*256 + (1:256);
  RFs(:, idx) = makeRF(Ms, randn(numel(zn), 256));
  RFr(:, idx) = makeRF(Mr, randn(numel(zn), 256));
end
if isfinite(snr_dB)
  zq = round(2*1.0/c*fs);  p0 = mean(mean(RFs(zq:zq+100, :).^2));
  RFs = RFs + sqrt(p0*10^(-snr_dB/10))*randn(size(RFs));
  RFr = RFr + sqrt(p0*10^(-snr_dB/10))*randn(size(RFr));
end
fprintf('Simulated %d + %d A-lines of %d scatterers in %.1f s\n', Nlines, Nlines, numel(zn), toc);

%% Gated power spectra (Hann gates, zero padded)
ng = round(2*Lg/c*fs);  T = ng/fs;  Np = 512;
w  = 0.5 - 0.5*cos(2*pi*((0:ng-1)' + 0.5)/ng);
fP = (0:Np/2)'*fs/Np;  kb = fP >= band(1) & fP <= band(2);  fb = fP(kb);  B = band(2) - band(1);
zc = (1.5:0.1:4.5)';                                           % gate centres (cm)
gateS = @(RF, z, cols) mean(abs(fft(RF(round(2*z/c*fs) - ng/2 + (1:ng), cols).*repmat(w, 1, numel(cols)), Np)).^2, 2);
Ss = zeros(Np, numel(zc));  Sr = Ss;
for i = 1:numel(zc)
  Ss(:, i) = gateS(RFs, zc(i), 1:Nlines);  Sr(:, i) = gateS(RFr, zc(i), 1:Nlines);
end
fprintf('Gate %.1f mm = %d samples = %.2f us; band %g-%g MHz; Hann; %d-point FFT\n', Lg*10, ng, T, band, Np);

%% 1. The factor 4: slope of the log ratio vs depth at each frequency (all lines, ROI 1.5-4.5 cm)
sel = zc >= 1.5 & zc <= 4.5;
RS  = log(Ss(kb, sel)./Sr(kb, sel));                            % ln power ratio, Np units
Zm  = [ones(sum(sel), 1) zc(sel)];
cf  = Zm \ RS';  slope_z = cf(2, :)';                           % d ln(Ss/Sr)/dz  per frequency (1/cm)
dA  = (alpha_s - alpha_r)/NP;                                   % true Delta alpha, Np/cm/MHz
fit1 = [ones(numel(fb), 1) fb] \ slope_z;
fprintf('\n1) d/dz ln(Ss/Sr) vs f: slope %.4f /cm/MHz; -4*dalpha = %.4f  -> ratio %.3f (expect 1)\n', ...
        fit1(2), -4*dA, fit1(2)/(-4*dA));
fprintf('   (power spectra, two-way: factor 4 = 2 (round trip) x 2 (power); amplitude spectra give 2)\n');

%% 2. Estimators on all 1024 lines
aRef = alpha_r - NP*fit1(2)/4;                                  % reference phantom method
RSn  = log(Ss(kb, sel));  cfn = Zm \ RSn';
fitn = [ones(numel(fb), 1) fb] \ cfn(2, :)';  aNoRef = -NP*fitn(2)/4;
i1 = find(abs(zc - 2) < 1e-9);  i2 = find(abs(zc - 4) < 1e-9);  dz = zc(i2) - zc(i1);
y2  = log(Ss(kb, i1)./Sr(kb, i1)) - log(Ss(kb, i2)./Sr(kb, i2));
p2  = [ones(numel(fb), 1) fb] \ y2;  aSLD = alpha_r + NP*p2(2)/(4*dz);
fprintf('\n2) alpha_s true %.3f dB/cm/MHz (reference %.2f)\n', alpha_s, alpha_r);
fprintf('   reference phantom method (regression 1.5-4.5 cm): %.3f dB/cm/MHz\n', aRef);
fprintf('   two-gate SLD, z = 2 and 4 cm, with reference : %.3f dB/cm/MHz\n', aSLD);
fprintf('   same regression WITHOUT reference            : %.3f dB/cm/MHz  <- diffraction bias\n', aNoRef);

%% 3. Variance: independent groups of K lines (reference averaged over all lines)
ng_ = Nlines/K;  aG = zeros(ng_, 1);  aR = aG;  aC = aG;
kc = fP > 0.5 & fP < 12;  sigp2 = sigf^2/2;                     % centroid band; power-spectrum variance (MHz^2)
fcr = (fP(kc)'*Sr(kc, :))./sum(Sr(kc, :));  pcr = Zm \ fcr(sel)'; % reference centroid slope vs depth (all lines)
SrA1 = Sr(kb, i1);  SrA2 = Sr(kb, i2);
zs = zc(sel);  j1 = find(abs(zs - 2) < 1e-9);  j2 = find(abs(zs - 4) < 1e-9);
for g = 1:ng_
  cols = (g-1)*K + (1:K);
  Sg = zeros(Np, numel(zs));
  for i = 1:numel(zs), Sg(:, i) = gateS(RFs, zs(i), cols); end   % each gate spectrum computed once
  y  = log(Sg(kb, j1)./SrA1) - log(Sg(kb, j2)./SrA2);           % two-gate SLD
  pg = [ones(numel(fb), 1) fb] \ y;  aG(g) = alpha_r + NP*pg(2)/(4*dz);
  cg = Zm \ log(Sg(kb, :)./Sr(kb, sel))';                       % reference method, all gates
  pr = [ones(numel(fb), 1) fb] \ cg(2, :)';  aR(g) = alpha_r - NP*pr(2)/4;
  Cg = (fP(kc)'*Sg(kc, :))./sum(Sg(kc, :));                      % centroid of every gate
  pcg = Zm \ Cg';  aC(g) = alpha_r - NP*(pcg(2) - pcr(2))/(4*sigp2);
end
psi1 = sum(1./(K:1e6).^2);                                      % trigamma psi1(K) without a toolbox
gw   = 1.75;                                                    % Hann-gate correlation factor (white-noise Monte Carlo)
sTh  = NP*sqrt(1.5*gw*psi1/(B^3*T))/dz;
fprintf('\n3) %d groups of K = %d lines:\n', ng_, K);
fprintf('   two-gate SLD : mean %.3f, std %.3f dB/cm/MHz;  formula sqrt(1.5 g_w psi1(K)/(B^3 T))/dz = %.3f\n', mean(aG), std(aG), sTh);
fprintf('   regression over %d overlapping gates, 1.5-4.5 cm: mean %.3f, std %.3f dB/cm/MHz\n', sum(sel), mean(aR), std(aR));
fprintf('   centroid, reference slope subtracted: mean %.3f, std %.3f dB/cm/MHz\n', mean(aC), std(aC));

%% 4. Centroid downshift on the same sample data (no reference)
fcen = (fP(kc)'*Ss(kc, :))./sum(Ss(kc, :));
pc = Zm \ fcen(sel)';  aCen = -NP*pc(2)/(4*sigp2);
aCenRef = alpha_r - NP*(pc(2) - pcr(2))/(4*sigp2);
fprintf('\n4) centroid slope %.4f MHz/cm; alpha = -slope/(4 sigma_p^2) = %.3f dB/cm/MHz (true %.2f)\n', pc(2), aCen, alpha_s);
fprintf('   expected slope without diffraction: -4 alpha sigma_p^2 = %.4f MHz/cm\n', -4*alpha_s/NP*sigp2);
fprintf('   with the reference centroid slope subtracted: alpha = %.3f dB/cm/MHz\n', aCenRef);

%% Figures
figure('Position', [100 100 1000 380]);
subplot(1, 3, 1);
zsel = zc(sel);  ic = find(ismember(round(10*zsel), [20 25 30 35 40]));
plot(fb, RS(:, ic), 'LineWidth', 1.2); grid on;
xlabel('frequency (MHz)'); ylabel('ln(S_s/S_r)  (Np)');
legend('z = 2.0 cm', '2.5', '3.0', '3.5', '4.0', 'Location', 'southwest');
title('Log spectral ratio vs depth');
subplot(1, 3, 2);
plot(fb, slope_z, 'b.', fb, fit1(1) + fit1(2)*fb, 'r-', fb, -4*dA*fb, 'k--', 'LineWidth', 1.2); grid on;
xlabel('frequency (MHz)'); ylabel('d ln(S_s/S_r)/dz  (1/cm)');
legend('per-frequency slope', 'linear fit', '-4\Delta\alpha f (truth)', 'Location', 'southwest');
title('Slope vs f gives -4\Delta\alpha');
subplot(1, 3, 3);
plot(zc, fcen, 'b-', zc, f0 - 4*alpha_s/NP*sigp2*zc, 'k--', 'LineWidth', 1.2); grid on;
xlabel('gate depth z (cm)'); ylabel('centroid (MHz)');
legend('measured (with focusing)', 'f_0 - 4\alpha\sigma_p^2 z', 'Location', 'southwest');
title('Centroid downshift');

What to look for: slope ratio to \(-4\Delta\alpha\) = 1.008; reference method 0.702 dB/cm/MHz (true 0.70) and 0.962 without the reference; two-gate SLD std 0.067 vs 0.069 predicted; regression std 0.030; raw centroid 1.13, reference-subtracted centroid 0.746 ± 0.030.

m27_bsc_esd_fit.mGaussian form factor fit in k², ESD sweep 20–300 µm showing the identifiable ka range, attenuation-compensation leakage · change band, K, z · runtime ≈ 6 s 
%% m27_bsc_esd_fit.m  -  backscatter coefficient shape, Gaussian form factor, effective scatterer diameter
% Purpose : speckle gates are simulated from a scatterer field whose power
%           spectrum is f^4 x F(f) with the Gaussian form factor
%           F = exp(-0.827 k^2 a^2) (a = effective radius, k = 2 pi f / c;
%           Insana et al. 1990), attenuated two-way to depth z, and seen through
%           a 5 MHz, 60 % system. A reference phantom with known BSC shape f^4
%           (Rayleigh scatterers) and known attenuation is scanned with the same
%           system. Then
%           (1) the reference phantom ratio gives BSC_s(f)/BSC_r(f) after
%               attenuation compensation exp(+4 (alpha_s - alpha_r) f z);
%           (2) ln(BSC_s / f^4) is linear in k^2 with slope -0.827 a^2, so a
%               linear fit returns the ESD = 2a;
%           (3) the ESD is identifiable only when ka over the band reaches ~0.5-1:
%               small scatterers give a nearly flat F and the fit fails;
%           (4) an attenuation-compensation error of 0.1 dB/cm/MHz at 3 cm
%               shifts the ESD by tens of micrometres (they share the k^2 / f slope).
% Change  : band, gate length Lg, lines per estimate K, depth z, alpha error.
% Expect  : ESD unbiased within 1 um and std 2-7 um once ka_max >= 1.4
%           (ESD >= 100 um at 7 MHz); ESD = 20 um: 36 % of fits give a slope of
%           the wrong sign and the rest average 33 um (biased, std 13 um); a
%           +0.1 dB/cm/MHz compensation error at z = 3 cm turns 100 um into 78 um.
clear; close all; rng(1);

%% Parameters
c   = 1540;                 % m/s  (= um/us, so k = 2 pi f / c in rad/um with f in MHz)
cc  = 0.154;                % cm/us, for gate lengths
fs  = 40;  f0 = 5;  bw = 0.6;
sigf = bw*f0/(2*sqrt(2*log(2)));
NP  = 20/log(10);           % 8.686 dB/Np
z   = 3.0;                  % gate depth (cm)
alpha_s = 0.5;  alpha_r = 0.5;   % dB/cm/MHz (equal here so attenuation errors are isolated)
Lg  = 0.4;  ng = round(2*Lg/cc*fs);  T = ng/fs;    % 4 mm gate
K   = 32;   Kref = 1024;    % sample lines per estimate, reference lines
band = [3 7];               % MHz
ntr = 200;                  % Monte Carlo estimates per ESD
Nf  = 1024;  ff = (0:Nf-1)'*fs/Nf;  ff(ff >= fs/2) = ff(ff >= fs/2) - fs;  fa = abs(ff);
w   = 0.5 - 0.5*cos(2*pi*((0:ng-1)' + 0.5)/ng);
Np  = 512;  fP = (0:Np/2)'*fs/Np;  kb = fP >= band(1) & fP <= band(2);  fb = fP(kb);  kk = 2*pi*fb/c;
G   = exp(-(fa - f0).^2/(2*sigf^2));
% amplitude filter for a gate: system x f^2 (Rayleigh) x sqrt(form factor) x two-way attenuation
Hamp = @(a_um, alphadB) G.*(fa/f0).^2.*exp(-0.827*(2*pi*fa/c*a_um).^2/2).*exp(-2*alphadB/NP*fa*z);
gateSpec = @(H, n) mean(abs(fft(real(ifft(fft(randn(Nf, n)).*repmat(H, 1, n)))...
                   (1:ng, :).*repmat(w, 1, n), Np)).^2, 2);

%% Reference phantom spectrum (known BSC shape f^4, known alpha_r), many lines
Sr = gateSpec(Hamp(0, alpha_r), Kref);  Sr = Sr(kb);
fprintf('Gate %.1f mm at z = %.1f cm, band %g-%g MHz, K = %d lines/estimate, reference %d lines\n', Lg*10, z, band, K, Kref);

%% 1-3. ESD sweep
ESDs = [20 40 60 80 100 150 200 300];          % um
fprintf('\n ESD(um)  ka range     mean ESD^  std    fails   ln-BSC slope n at f0\n');
X = [ones(numel(kk), 1) kk.^2];
keep = zeros(numel(fb), 3);
for e = 1:numel(ESDs)
  a = ESDs(e)/2;  est = nan(ntr, 1);
  for t = 1:ntr
    Ss = gateSpec(Hamp(a, alpha_s), K);
    R  = Ss(kb)./Sr.*exp(4*(alpha_s - alpha_r)/NP*fb*z);   % = BSC_s/BSC_r = F(f) x const
    p  = X \ log(R);                                       % ln R = b0 - 0.827 a^2 k^2
    a2 = -p(2)/0.827;
    if a2 > 0, est(t) = 2*sqrt(a2); end
    if t == 1 && any(ESDs(e) == [20 100 300]), keep(:, find(ESDs(e) == [20 100 300])) = log(R) - p(1); end
  end
  n_f0 = 4 - 2*0.827*(2*pi*f0/c*a)^2;                     % local power-law exponent of f^4 F at f0
  fprintf('  %4d    %.2f-%.2f   %7.1f  %6.1f   %5.3f   %6.2f\n', ESDs(e), kk(1)*a, kk(end)*a, ...
          mean(est(~isnan(est))), std(est(~isnan(est))), mean(isnan(est)), n_f0);
end
fprintf('BSC = f^4 exp(-0.827 k^2 a^2) peaks at ka = sqrt(2/0.827) = %.3f; local exponent n = 4 - 1.654 (ka)^2\n', sqrt(2/0.827));

%% 4. Attenuation-compensation error leaks into the ESD (noise-free spectra)
fprintf('\nCompensation error (assumed alpha_s - true), z = %.1f cm, noise-free expected spectra:\n', z);
Hn = @(a) (Hamp(a, alpha_s)).^2;                        % expected power (no window smoothing)
for a_true = [50 100]
  for dA = [-0.1 0.05 0.1 0.2]
    Pn = Hn(a_true);  Pr = Hamp(0, alpha_r).^2;
    Pn = interp1(ff(1:Nf/2), Pn(1:Nf/2), fb);  Pr = interp1(ff(1:Nf/2), Pr(1:Nf/2), fb);
    R  = Pn./Pr.*exp(4*(alpha_s + dA - alpha_r)/NP*fb*z);
    p  = X \ log(R);  a2 = -p(2)/0.827;
    if a2 > 0, s = sprintf('%.1f um', 2*sqrt(a2)); else, s = 'no real solution'; end
    fprintf('  true ESD %3d um, error %+5.2f dB/cm/MHz -> ESD^ = %s\n', 2*a_true, dA, s);
  end
end

%% Figures
figure('Position', [100 100 950 380]);
subplot(1, 2, 1);
kaa = linspace(0, 4, 400);
plot(kaa, 10*log10(kaa.^4.*exp(-0.827*kaa.^2)/max(kaa.^4.*exp(-0.827*kaa.^2))), 'k', 'LineWidth', 1.4); hold on;
plot(kaa, 10*log10(exp(-0.827*kaa.^2)), 'b--', 'LineWidth', 1.2);
for a = [10 50 150], r = 2*pi*band/c*a; plot(r, [-29 -29], 'LineWidth', 4); end
ylim([-30 2]); grid on; xlabel('ka'); ylabel('dB');
legend('(ka)^4 F(ka), normalized', 'F = exp(-0.827 (ka)^2)', '3-7 MHz, ESD 20 um', 'ESD 100 um', 'ESD 300 um', 'Location', 'southeast');
title('Gaussian form factor and the band in ka');
subplot(1, 2, 2);
plot(kk.^2*1e3, keep, '.', 'MarkerSize', 8); hold on;
for a = [10 50 150], plot(kk.^2*1e3, -0.827*a^2*kk.^2, 'k-'); end
grid on; xlabel('k^2  (10^{-3} rad^2/um^2)'); ylabel('ln(BSC_s/BSC_r) - intercept  (Np)');
legend('ESD 20 um', 'ESD 100 um', 'ESD 300 um', 'Location', 'southwest');
title('One estimate (K lines): slope = -0.827 a^2');

What to look for: ESD 20 µm: 36 % failed fits, survivors 33 ± 13 µm; 100 µm: 99 ± 7 µm; 300 µm: 296 ± 2 µm; compensation error +0.10 dB/cm/MHz at 3 cm turns 100 µm into 77.6 µm and 200 µm into 189.8 µm.

27.14Self-check

Q1A single Hann-windowed periodogram of a 256-sample speckle gate: what are the mean and standard deviation of \(10\log_{10}(P/S)\) in band?

−2.51 dB and 5.57 dB (27.5), independent of \(N\) and of the window (the window changes the bias, not the \(\chi^2_2\) distribution).

Q2Welch with N = 1024, M = 128, Hann, 50 % overlap: how many segments, how many effective looks, and how wide is the 95 % interval?

\(N_b = 15\); \(K_{\mathrm{eff}} = 15/(1 + 2\cdot(14/15)\cdot(1/6)^2) = 14.3\); \(\nu = 28.5\); interval −1.99 to +2.60 dB around \(\hat S\), a 4.6 dB span.

Q3\(\ln S(z_2) - \ln S(z_1)\) has slope −0.46 per MHz for gates 2 cm apart (power spectra). What is \(\alpha_0\)? What if these were amplitude spectra?

Power: \(\alpha_0 = 0.46/(4\cdot 2)\) = 0.0575 Np/cm/MHz = 0.50 dB/cm/MHz. Amplitude: \(\alpha_0 = 0.46/(2\cdot 2)\) = 0.115 Np/cm/MHz = 1.0 dB/cm/MHz (the same data would have shown slope −0.23 in amplitude).

Q4N = 128 at 40 MHz with NW = 4: how many tapers, what resolution, what variance reduction?

\(K = 2NW - 1 = 7\) tapers, resolution \(2W = 2NW f_s/N\) = 2.5 MHz, variance about 1/7 of a single periodogram (0.17 measured on speckle because the spectrum is not flat over 2.5 MHz).

Q5A 5 MHz, 60 % Gaussian pulse shows a centroid slope of −0.187 MHz/cm. What attenuation does it imply, and what must you assume?

\(\sigma_p^2 = \sigma_f^2/2 = 1.274^2/2\) = 0.811 MHz²; \(\alpha_0 = 0.187/(4\cdot 0.811)\) = 0.0576 Np/cm/MHz = 0.50 dB/cm/MHz. Assumes a Gaussian power spectrum unchanged in shape by BSC, diffraction and noise, and \(\alpha \propto f\).

Q6Can you estimate an ESD of 20 µm with a 3–7 MHz probe? Of 150 µm?

20 µm: \(ka\) = 0.12–0.29, local exponent ≈ 3.9: no, the form factor is flat within the noise. 150 µm: \(ka\) = 0.92–2.14, exponent at 5 MHz ≈ 0.1: yes, well identified (149 ± 4 µm in the script).

Q7You double the ROI length of an attenuation estimate from 1.5 to 3 cm and halve the number of lines. Does precision improve?

\(\sigma \propto K^{-1/2}L^{-3/2}\): factor \(\sqrt{2}\) worse from \(K\), \(2^{1.5}\) = 2.83 better from \(L\), net 2.0× better.

27.15Figures from the literature

Published figures that put measured data behind the estimators and numbers in this module. Citations and figure numbers were checked against the papers; images are not reproduced, follow the links.

The two calibration routes for removing the system and the beam

M. L. Oelze and J. Mamou, "Review of quantitative ultrasound: envelope statistics and backscatter coefficient imaging and contributions to diagnostic ultrasound," IEEE TUFFC 63(2):336–351, 2016, Fig. 1. doi:10.1109/TUFFC.2015.2513958

Fig. 1 contrasts the planar-reflector calibration (a reference echo from a flat reflector, then an analytic diffraction correction) with the reference-phantom calibration (a phantom of known \(\alpha_r\) and \(\mathrm{BSC}_r\) scanned with identical settings). Both exist to remove \(G\) and \(D\) from (27.16); the reference phantom of Section 27.5 also removes depth-dependent focusing and TGC without modelling them, which is why it dominates array-based QUS.

ESD images of real tumours, and a calibrated BSC that tracks liver fat

Same review, Figs. 2 and 8.

Fig. 2 overlays ESD estimates on B-mode images of rat fibroadenomas, mouse carcinomas and mouse sarcomas: the Gaussian-form-factor fit of (27.24) applied per ROI. Fig. 8 shows measured BSC versus frequency for liver samples after 0, 3 and 6 weeks of a fatty diet. This is the absolute, calibrated quantity (27.23) delivers once \(G\), \(D\) and the attenuation path are removed; ask of any such plot how many independent lines and which attenuation compensation went into each curve.

Two laboratories, the same phantoms: how well do BSC and attenuation agree?

J. J. Anderson, M.-T. Herd, M. R. King, A. Haak, Z. T. Hafez, J. Song, M. L. Oelze, E. L. Madsen, J. A. Zagzebski, W. D. O'Brien Jr. and T. J. Hall, "Interlaboratory comparison of backscatter coefficient estimates for tissue-mimicking phantoms," Ultrasonic Imaging 32(1):48–64, 2010, Figs. 6 and 12. doi:10.1177/016173461003200104

Fig. 6: insertion-loss attenuation of the same phantoms measured independently in two laboratories (agreement within about 0.01 dB/cm/MHz mean-square difference). Fig. 12: BSC versus frequency for a phantom of 41 µm glass spheres from both laboratories against the Faran-theory prediction; the paper reports average fractional errors of a few percent there, about 8–14 % averaged over all phantoms, and an estimate uncertainty near 30 %. That 30 % is about 1.1 dB: compare it with the \(\chi^2\) interval of (27.8) for the number of looks you average.

Four clinical array systems, one phantom: reference-phantom attenuation and ESD in practice

K. Nam, I. M. Rosado-Mendez, L. A. Wirtzfeld, A. D. Pawlicki, V. Kumar, E. L. Madsen, G. Ghoshal, R. J. Lavarello, M. L. Oelze, T. A. Bigelow, J. A. Zagzebski, W. D. O'Brien Jr. and T. J. Hall, "Ultrasonic attenuation and backscatter coefficient estimates of rodent-tumor-mimicking structures: comparison of results among clinical scanners," Ultrasonic Imaging 33(4):233–250, 2011, Figs. 4–6 and Tables 5 and 7. doi:10.1177/016173461103300403

Figs. 4 and 5: BSC versus frequency of two glass-bead inclusions measured with laboratory single-element transducers and with Siemens, Ultrasonix, VisualSonics and Zonare array systems, against Faran predictions. Fig. 6: per-system percentage differences. Table 5: attenuation from the sample-to-reference spectral ratio versus depth (the method of (27.19)) agreed with the laboratory values to 11 % on average and 29 % at worst. Table 7: ESD scatter between systems is larger in the low band than in the high band (e.g. ±12 µm vs ±1 µm for the 35 µm beads), the identifiability rule of Section 27.7 seen in measured data.

27.16Go deeper

  • D. B. Percival & A. T. Walden, Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press, 1993. Periodogram statistics, Welch, DPSS and the tridiagonal method.
  • D. J. Thomson, "Spectrum estimation and harmonic analysis," Proceedings of the IEEE, 1982.
  • D. Slepian, "Prolate spheroidal wave functions, Fourier analysis, and uncertainty V: the discrete case," Bell System Technical Journal, 1978.
  • P. D. Welch, "The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms," IEEE Trans. Audio and Electroacoustics, 1967.
  • S. M. Kay, Modern Spectral Estimation: Theory and Application, Prentice Hall, 1988. AR, Burg, Yule–Walker, order selection.
  • H. Akaike, "A new look at the statistical model identification," IEEE Trans. Automatic Control, 1974.
  • L. X. Yao, J. A. Zagzebski & E. L. Madsen, "Backscatter coefficient measurements using a reference phantom to extract depth-dependent instrumentation factors," Ultrasonic Imaging 12(1):58–70, 1990.
  • M. F. Insana, R. F. Wagner, D. G. Brown & T. J. Hall, "Describing small-scale structure in random media using pulse-echo ultrasound," J. Acoust. Soc. Am. 87(1):179–192, 1990.
  • F. L. Lizzi, M. Greenebaum, E. J. Feleppa, M. Elbaum & D. J. Coleman, "Theoretical framework for spectrum analysis in ultrasonic tissue characterization," J. Acoust. Soc. Am., 1983.
  • M. L. Oelze & J. Mamou, "Review of quantitative ultrasound: envelope statistics and backscatter coefficient imaging and contributions to diagnostic ultrasound," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63(2):336–351, 2016.
  • Y. Labyed & T. A. Bigelow, "A theoretical comparison of attenuation measurement techniques from backscattered ultrasound echoes," J. Acoust. Soc. Am., 2011.
  • A. L. Coila & R. Lavarello, "Regularized spectral log difference technique for ultrasonic attenuation imaging," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 65(3), 2018.