Part IV · Module 24

Fourier-domain beamforming: f-k migration and k-space imaging

Module 23 read every transmit scheme as a way of filling k-space. This module turns that picture into an imaging algorithm: transform the channel data once, move every spectral sample to the object wavenumber it measures, transform back. Stolt's f-k migration and its plane-wave versions are the result: the physics of DAS at FFT cost, with every approximation visible as a mapping error you can plot and bound.

You will be able to

Read channel data in \((k_x, f)\): the propagating cone, the evanescent region, spatial aliasing for plane-wave and zero-offset data
Derive Stolt's map for zero-offset data from the exploding-reflector model, including the \(c/2\) and the Jacobian, and know its one homogeneous-medium error (interpolation)
Derive Garcia's plane-wave ERM fit \((\hat{v}, \beta, \gamma)\), bound its error, and write the exact transmit-plus-receive wavevector map
Predict the k-space support and PSF of a plane-wave set and compound in k-space
Defend f-k against DAS on accuracy (F#, aliasing, interpolation, gain), cost and the medium model, and extend it to layers by phase shift

Assumes

Module 04 · Sampling, Fourier & LTI systems: 2D FFT, k-space of an image, beam = FT of aperture
Module 08 · Analytic signal & IQ demodulation: IQ data and the analytic signal; DAS with steered plane waves
Module 23 · Beam patterns, k-space and compounding: coarray, plane-wave compounding, the \(\lambda/\Delta(\sin\theta)\) replica rule

Exact formulaclosed form plotted as is (dispersion cone, Stolt and Ewald maps, ERM parameters)SimulationJS in the page: point and speckle channel data synthesized from the delay law, FFT-based f-k migration and DAS on the same dataApprox. modelsimplified physics (narrowband PSF, flop counts, straight-ray layer shifts, Monte Carlo coefficients); assumptions stated under the badgeSchematicdrawing for intuition; the k-space drawings are to scale at 5 MHz but not computed from dataAll simulations are our own code (JS here, Octave/MATLAB in 24.11); no Field II or k-Wave output is shown. c = 1540 m/s throughout.

24.1The angular spectrum: channel data as a sum of plane waves

In a homogeneous medium every solution of the wave equation is a superposition of plane waves \(e^{\,j2\pi(k_x x + k_z z + f t)}\) (sign of \(f t\) chosen to match the transform in (24.1), so positive f is the analytic band). Wavenumbers are in cycles/mm and frequencies in MHz throughout, with \(c = 1.54\) mm/µs. The wave equation ties the three together, \(k_x^2 + k_z^2 = (f/c)^2\), so the array, which measures only x and t, determines \(k_z\) up to a sign:

\[ \begin{aligned} D(k_x, f) &= \iint d(x, t)\, e^{-j2\pi(k_x x + f t)}\, dx\, dt \\ k_z &= \pm\sqrt{(f/c)^{2} - k_x^{2}} \\ k_x &= (f/c)\sin\alpha \end{aligned} \]

(24.1)

\(d(x, t)\): channel data: element position x (mm), fast time t (µs)
\(D(k_x, f)\): its 2D spectrum (MATLAB: fft over t, then over x)
\(k_x, k_z\): lateral and axial wavenumbers, cycles/mm
\(f\): temporal frequency, MHz
\(c\): speed of sound, 1.54 mm/µs
\(\alpha\): arrival angle to the array normal
\(\iint\): double integral over element position and time
\(j\): imaginary unit
\(x,\ t\): as in the definition of \(d(x, t)\)

\(k_z\) is real for \(\lvert k_x\rvert \le \lvert f\rvert/c\) (propagating) and imaginary otherwise (evanescent). A sample at \((k_x, f)\) is a plane wave arriving from angle α; the cone edge \(\lvert k_x\rvert = f/c\) is a grazing arrival. At 5 MHz, \(f/c = 3.25\) cycles/mm.

Where the data live. In plane-wave imaging the receive leg is a one-way wave from each scatterer to the array, so the data spectrum fills the cone \(\lvert k_x\rvert \le f/c\). Zero-offset (monostatic) data have delay \(2r/c\) and slope \(2\sin\alpha/c\) across the array, so they fill \(\lvert k_x\rvert \le 2f/c\). Outside the cone the field would be evanescent:

\[ \begin{gathered} \lvert p(z)\rvert \propto e^{-2\pi z\sqrt{k_x^{2} - (f/c)^{2}}} \\ 5\ \mathrm{MHz}:\ -25,\ -81,\ -198\ \mathrm{dB/mm} \\ \text{at}\ k_x = 1.01,\ 1.1,\ 1.5 \times f/c \end{gathered} \]

(24.2)

\(\lvert p(z)\rvert\): amplitude of an evanescent component at distance z from the source plane
\(z\): distance from the source plane, mm
\(k_x\): lateral wavenumber, cycles/mm
\(f/c\): wavenumber of a plane wave at frequency \(f\) (radius of the propagating cone), cycles/mm
\(f,\ c\): as in (24.1)
\(\propto\): proportional to
\(\mathrm{dB/mm}\): 20 log\(_{10}\) amplitude loss per mm of distance

Scatterers sit millimetres from the array, so their evanescent components arrive more than 100 dB down. What appears outside the cone in practice is leakage from truncating the aperture: −22 dB for a target at 5 mm with a rectangular aperture, −29 to −53 dB with a 25 % Tukey taper (m24_kspace_evanescent.m).

The phase-shift propagator. Moving the recorded field from z to \(z + \Delta z\) in a homogeneous slab multiplies each component by a pure phase. This is the angular spectrum method, exact for propagating waves:

\[ \begin{aligned} P(k_x, f; z + \Delta z) &= P(k_x, f; z) \\ &\quad \times e^{\,j2\pi k_z \Delta z} \\ k_z &= \sqrt{(f/c)^{2} - k_x^{2}} \\ &\quad (\lvert k_x\rvert \le f/c\ \text{only}) \end{aligned} \]

(24.3)

\(P(k_x, f; z)\): angular spectrum of the field on the plane at depth z
\(\Delta z\): propagation step, mm
\(k_x,\ f,\ k_z\): as in (24.1)
\(z\): depth, mm
\(c\): speed of sound, mm/µs
\(j\): imaginary unit

Back-propagating evanescent components would multiply them by \(e^{+2\pi\lvert k_z\rvert\Delta z}\): 81 dB per millimetre at \(1.1\,f/c\), applied to noise. Every Fourier imaging method masks them.

Spatial sampling. Elements at pitch p sample x, so \(k_x\) is periodic with period \(1/p\) and the Nyquist limit is \(1/(2p)\). The plane-wave cone stays inside it at every angle only if \(p \le \lambda_{\min}/2\). A 0.30 mm pitch is about λ at 5 MHz: the cone meets Nyquist at \(f = c/(2p) = 2.57\) MHz, and above it every arrival steeper than \(\alpha_{\mathrm{alias}} = \arcsin\!\big(c/(2pf)\big)\) folds (30.9° at 5 MHz, 20.0° at 7.5 MHz). For a 0° plane wave the share of channel energy beyond \(1/(2 \times 0.30\ \mathrm{mm})\) is 82 % for a target at 5 mm, 65 % at 10 mm, 34 % at 20 mm and 10 % at 30 mm.

Demo 24.1

Channel data and its k_x–f spectrum: cone, Nyquist and folding

SimulationIQ channel data synthesized in JS from the plane-wave or zero-offset delay law (Gaussian pulse, point targets), 2D FFT in the page; cone and Nyquist lines are exact.

cone \(\lvert k_x\rvert = f/c\) (×2 for zero-offset)Nyquist \(\pm 1/(2p)\)cone folded by \(\pm 1/p\)

Data

Pitch 0.30 mm

Target depth 12 mm

Steering θ 0°

steering compensation

Try: at pitch 0.075 mm everything sits inside the cone and inside Nyquist. Step the pitch up: at 0.30 mm the red folded copies enter the band above 2.6 MHz. Move the target to 4 mm (steep arrivals, more folding) and to 30 mm (shallow arrivals, little). Switch to zero-offset: the cone doubles and even 0.15 mm aliases. Steer to 15° and tick compensation: the cone shears by \(f\sin\theta/c\), the operation every plane-wave f-k method starts with (Section 24.3). 19.2 mm aperture, 5 MHz, 60 % band, IQ at 20 MHz, two point targets.

Pitfall

The receive cone is not the image's k-space. The image spectrum of a plane-wave acquisition is the receive wavevector plus the transmit wavevector (Section 24.4), which doubles the axial wavenumber and shifts the lateral one by \((f/c)\sin\theta\). Reading resolution off the raw data spectrum underestimates the axial extent by a factor of two.

24.2Exploding reflectors and Stolt's change of variables

The exploding reflector. In zero-offset data each element fires alone and records its own echo. A point at \((x_0, z_0)\) arrives at

\[ \begin{aligned} \tau(x) &= \frac{2\sqrt{(x - x_0)^{2} + z_0^{2}}}{c} \\ &= \frac{\sqrt{(x - x_0)^{2} + z_0^{2}}}{c/2} \end{aligned} \]

(24.4)

\((x_0, z_0)\): scatterer position, mm
\(x\): position of the element that fires and receives, mm
\(\tau(x)\): zero-offset (monostatic) echo delay, µs
\(c\): speed of sound, mm/µs
\(c/2\): effective speed of the exploding-reflector model

This is the one-way arrival time of a source at the reflector that fires at t = 0 into a medium of speed \(v = c/2\). Every reflector "explodes" at t = 0, the recorded data are a one-way wave field, and imaging means propagating it back to t = 0 (Loewenthal et al., 1976; Claerbout, 1985).

Downward continuation by phase shift gives the field at depth z; the imaging condition t = 0 sums it over frequency. Gazdag's phase-shift migration (1978) evaluates this depth by depth:

\[ \begin{aligned} r(x, z) &= \iint D(k_x, f)\, e^{\,j2\pi(k_x x + k_z z)}\, df\, dk_x \\ k_z &= \operatorname{sign}(f)\sqrt{(f/v)^{2} - k_x^{2}} \\ v &= c/2 \end{aligned} \]

(24.5)

\(r(x, z)\): complex image (reflectivity estimate)
\(v\): migration speed, \(c/2\) for zero-offset data
\(D(k_x, f)\): data spectrum, as in (24.1)
\(k_x,\ k_z,\ f\): as in (24.1)
\(\operatorname{sign}(f)\): +1 for positive and \(-1\) for negative frequencies
\(\iint\): double integral over \(f\) and \(k_x\)
\(c\): speed of sound, mm/µs

Sign check: a vertical upgoing wave seen at depth z arrives \(z/v\) earlier, a factor \(e^{+j2\pi f z/v}\) with this transform convention, so \(k_z = +f/v\) for \(f > 0\).

Stolt's step (Stolt, 1978) removes the depth loop: change the integration variable from f to \(k_z\). The integral becomes an inverse 2D Fourier transform of a resampled, reweighted copy of the data spectrum:

\[ \begin{aligned} f &= v\,\operatorname{sign}(k_z)\sqrt{k_x^{2} + k_z^{2}} \\ df &= \frac{v\,\lvert k_z\rvert}{\sqrt{k_x^{2} + k_z^{2}}}\, dk_z \\ R(k_x, k_z) &= D\big(k_x,\ v\sqrt{k_x^{2} + k_z^{2}}\big) \\ &\quad \times \frac{v\,k_z}{\sqrt{k_x^{2} + k_z^{2}}} \\ r &= \mathrm{IFFT}_2\left[R\right] \end{aligned} \]

(24.6)

\(R(k_x, k_z)\): image spectrum on the output wavenumber grid
\(v k_z/\sqrt{k_x^2 + k_z^2}\): Jacobian \(\lvert\partial f/\partial k_z\rvert = v\cos\varphi\), φ the propagation angle
\(\mathrm{IFFT}_2\): inverse 2D FFT over \((k_x, k_z)\)
\(f,\ k_x,\ k_z\): as in (24.1) and (24.5)
\(D\): data spectrum, as in (24.1)
\(df\): frequency element mapped to \(dk_z\) by the change of variables (Stolt mapping)
\(v\): migration speed, as in (24.5)

The weight is the Jacobian: the cosine of the propagation angle of the plane wave. Since \(f \ge v\lvert k_x\rvert\) on the output grid, the evanescent region is never read. Only \(f > 0\) is used, so the image is analytic and \(\lvert r\rvert\) is the envelope.

FFT in t and in x, zero-padded (Section 24.5 says how much).
For each \(k_x\) column and each output \(k_z\), compute \(f = v\sqrt{k_x^2 + k_z^2}\) and interpolate \(D(k_x, \cdot)\) there: one 1D interpolation per column.
Multiply by the Jacobian.
Inverse 2D FFT.

What it does to a sample. A sample at \((k_x, f)\) moves to \(k_z = \sqrt{(f/v)^2 - k_x^2}\): straight down in \(k_z\), along its \(k_x\) column (Figure 24.A). A flat reflector (\(k_x = 0\)) is only relabelled, \(k_z = 2f/c\), which is time-to-depth conversion. Steep components (large \(k_x\)) move most; this is the dip correction that collapses hyperbolas.

Schematic Figure 24.A · Stolt's map in k-space: each sample is first placed at k_z = 2f/c and then moved straight down its k_x column onto the circle of radius 2f/c (drawn to scale at 5 MHz, not computed from data)

Interpolation is the only error in a homogeneous medium. Linear interpolation of \(e^{-j2\pi f t}\) between bins Δf apart at fractional position w is off by \(w(1 - w)\varphi^2/2\) with \(\varphi = 2\pi\,\Delta f\, t\). The mean part is a smooth gain loss (\(\varphi^2/12\)); the part that varies from one \(k_x\) column to the next (rms \(0.037\,\varphi^2\), verified by Monte Carlo) spreads each echo across x at its own depth. Two cures, both cheap: zero-pad the record before the FFT (smaller Δf) and centre the time origin (multiply by \(e^{\,j2\pi f t_c}\) before, \(e^{-j2\pi F t_c}\) after, so t becomes \(t - t_c\)).

\[ \begin{gathered} \text{smear} \approx 20\log_{10}\!\left(0.037\,\varphi^{2}\right)\ \mathrm{dB\ per\ echo} \\ \varphi = 2\pi\,\Delta f\,\lvert t - t_c\rvert \\ -30\ \mathrm{dB\ needs}\ \varphi \le 0.92 \end{gathered} \]

(24.7)

\(\Delta f\): FFT bin spacing, \(f_s/N_{\mathrm{FFT}}\)
\(t_c\): time origin after centring, µs
\(\varphi\): phase step of the echo between adjacent bins, rad
\(\text{smear}\): level of the interpolation artifact relative to the echo, dB
\(\log_{10}\): base-10 logarithm
\(t\): echo time, µs
\(\lvert\cdot\rvert\): absolute value

m24_stolt_planewave.m, 62.5 µs record: a 2048-point FFT without centring (\(\varphi \approx 3\) at the speckle) leaves clutter beside a speckle block at −10.6 dB instead of the DAS −18.7 dB; 4096 points give −18.2 dB, and centring adds headroom.

Velocity and sampling. Migrating zero-offset data with \(v = c\) instead of \(c/2\) halves every \(k_z\): a target at 10 mm reappears at 19.94 mm, 21 to 24 dB below its focused peak (m24_fk_migration_zero_offset.m). Zero-offset data fill \(\lvert k_x\rvert \le 2f/c\), so the alias-free pitch is \(\lambda_{\min}/4\) (0.077 mm at 5 MHz, 0.048 mm at 8 MHz), or \(\lambda/(4\sin\alpha)\) once a receive F-number limits the angle.

Demo 24.2

Stolt mapping, step by step (zero-offset data, exploding reflector)

SimulationZero-offset IQ data of point targets synthesized in JS (256 elements, 0.075 mm), Stolt map with the selected interpolation, padding and Jacobian, inverse FFT in the page.

output row \(k_z\) probe: where it samples Dtrue target positions

Step 2 · map

Probe output row k_z 6.5 /mm

Migration speed v/c 0.50

Interpolation along f

FFT length / record

Jacobian weight centre the time origin

Try: slide the \(k_z\) probe: each output row reads the data spectrum along a hyperbola \(f = v\sqrt{k_x^2 + k_z^2}\) that touches \(k_x = 0\) at \(f = v k_z\). Set v/c to 1.00: targets land at twice their depth and smear into arcs. Nearest-bin interpolation with no padding fills the image with lateral streaks at the target depths; linear, 2×, centred removes them. Turning the Jacobian off barely changes the image: it is an amplitude weight, not a focusing term. 256 elements at 0.075 mm (λ/4 at 5 MHz), fan filter F# 1.

Intuition

DAS sums along a hyperbola in (x, t) for every pixel. Stolt does the same bookkeeping once for all pixels: in the Fourier domain a hyperbola family becomes a family of curves \(f(k_x, k_z)\), and summing along all of them at once is a resampling followed by an inverse FFT.

24.3Plane-wave transmits: steering compensation and Garcia's ERM fit

For a plane wave steered by θ, with the wavefront through the array centre at t = 0, a point at \((x_0, z_0)\) reaches element x at

\[ \begin{aligned} c\,\tau(x) &= x_0\sin\theta + z_0\cos\theta \\ &\quad + \sqrt{(x - x_0)^{2} + z_0^{2}} \end{aligned} \]

(24.8)

\(\theta\): steering angle of the transmitted plane wave
\(x_0\sin\theta + z_0\cos\theta\): transmit path to the scatterer, mm
\(\sqrt{(x - x_0)^2 + z_0^2}\): receive path to element x, mm
\(\tau(x)\): plane-wave echo delay at element x, µs
\(c\): speed of sound, mm/µs
\(x,\ x_0,\ z_0\): element and scatterer positions, as in (24.4), mm

The transmit term depends on \(x_0\), the receive term on \(x - x_0\). A migration needs a curve that depends only on \(x - x_0\).

Steering compensation. Shift every channel by \(-x\sin\theta/c\), i.e. multiply \(D(x, f)\) by \(e^{\,j2\pi f x\sin\theta/c}\). With \(u = x - x_0\) the curve becomes shift-invariant:

\[ \begin{gathered} c\,\tau'(u) = z_0\cos\theta - u\sin\theta + \sqrt{u^{2} + z_0^{2}} \\ \Longleftrightarrow\quad k_x \to k_x + \frac{f\sin\theta}{c} \end{gathered} \]

(24.9)

\(\tau^{\prime}(u)\): steering-compensated echo delay, µs
\(u\): \(x - x_0\), lateral offset between element and scatterer, mm
\(z_0,\ \theta,\ c\): as in (24.8)
\(\Longleftrightarrow\): equivalent to
\(k_x \to k_x + f\sin\theta/c\): shear of the lateral wavenumber axis (cycles/mm)
\(f\): temporal frequency, MHz

In the spectrum this is a shear: the receive cone tilts by \(\sin\theta\) (Demo 24.1). On a sampled array the sheared spectrum can exceed Nyquist, so upsample laterally first (the MATLAB code zero-pads \(k_x\) by 2×) or the image wraps.

\(\tau'\) is not a hyperbola, so Stolt cannot be applied directly. Garcia et al. (2013) fit an exploding-reflector hyperbola \(\hat{\tau}(u) = \sqrt{(u - \hat{u})^2 + \hat{z}^2}/\hat{v}\) that matches \(\tau'\) in value, slope and curvature at u = 0, the element above the scatterer:

\[ \begin{aligned} \frac{R}{\hat{v}} &= \frac{z_0(1 + \cos\theta)}{c} &&\text{value} \\ \frac{\hat{u}}{R\hat{v}} &= \frac{\sin\theta}{c} &&\text{slope} \\ \frac{\hat{z}^{2}}{R^{3}\hat{v}} &= \frac{1}{z_0 c} &&\text{curv.} \\ R^{2} &= \hat{u}^{2} + \hat{z}^{2} \\ \Rightarrow\ R^{2} &= \frac{z_0^{2}(1 + \cos\theta)^{2}}{1 + \cos\theta + \sin^{2}\theta} \end{aligned} \]

(24.10)

\(\hat{v}\): equivalent exploding-reflector speed, mm/µs
\((\hat{u}, \hat{z})\): apex of the fitted hyperbola relative to the scatterer column, mm
\(R\): distance from the element above the scatterer to the fitted apex, mm
\(z_0,\ \theta,\ c\): as in (24.8)

\[ \begin{aligned} \hat{v} &= \frac{c}{\sqrt{1 + \cos\theta + \sin^{2}\theta}} \\ \beta &= \frac{(1 + \cos\theta)^{3/2}}{1 + \cos\theta + \sin^{2}\theta}, \quad \hat{z} = \beta z_0 \\ \gamma &= \frac{\sin\theta}{2 - \cos\theta}, \quad \hat{u} = \gamma z_0 \end{aligned} \]

(24.11)

\(\hat{v}\): ERM migration speed; \(c/\sqrt{2}\) at θ = 0
\(\beta\): depth stretch of the ERM image, \(\sqrt{2}\) at θ = 0
\(\gamma\): lateral shear of the ERM image per unit depth, 0 at θ = 0
\(\theta\): steering angle
\(c\): speed of sound, mm/µs
\(z_0\): scatterer depth, mm

θ = 0: \(\hat{v} = c/\sqrt{2}\), \(\beta = \sqrt{2}\), \(\gamma = 0\). θ = 10°: \(\hat{v} = 0.7045c\), β = 1.388, γ = 0.171. θ = 20°: 0.6973c, 1.314, 0.323. The identity \((1 + \cos\theta)/(1 + \cos\theta + \sin^2\theta) = 1/(2 - \cos\theta)\) gives γ its compact form.

Back to true coordinates. Migration with \(\hat{v}\) puts the point at \((x_0 + \gamma z_0, \beta z_0)\). The image we want is \(I(x, z) = \hat{I}(x + \gamma z, \beta z)\); in k-space that is the substitution \(\hat{k}_z = (K_z - \gamma K_x)/\beta\), so the whole method is still one interpolation along f per column:

\[ \begin{aligned} f &= \hat{v}\sqrt{K_x^{2} + \left(\frac{K_z - \gamma K_x}{\beta}\right)^{2}} \\ \left\lvert\frac{\partial f}{\partial K_z}\right\rvert &= \frac{\hat{v}\,\hat{k}_z}{\beta\sqrt{K_x^{2} + \hat{k}_z^{2}}} \end{aligned} \]

(24.12)

\((K_x, K_z)\): object (image) wavenumbers on the output grid, cycles/mm
\(\hat{k}_z\): axial wavenumber of the ERM image, \((K_z - \gamma K_x)/\beta\)
\(f\): temporal frequency read from the data for the output sample, MHz
\(\hat{v},\ \beta,\ \gamma\): as in (24.11)
\(\partial f/\partial K_z\): Jacobian of the change of variables, mm/µs
\(\lvert\cdot\rvert\): absolute value

The price of the fit. It is exact to second order at u = 0. At θ = 0 the leading error is a path excess of \(z_0\tan^4\alpha/16\) at receive angle α (exact: 21 µm at 18.4° and 99 µm at 26.6° for \(z_0 = 30\) mm; the formula gives 23 and 117 µm). Keeping it below λ/4, a quarter-cycle of phase error, bounds the usable receive aperture:

\[ \begin{gathered} F\# \ge \left(\frac{z}{64\lambda}\right)^{1/4} \\ 5\ \mathrm{MHz}:\ 1.00,\ 1.11,\ 1.19,\ 1.32 \\ \text{at}\ z = 20,\ 30,\ 40,\ 60\ \mathrm{mm} \end{gathered} \]

(24.13)

\(F\#\): receive f-number, depth over aperture width, dimensionless
\(z\): imaging depth, mm
\(\lambda\): wavelength \(c/f\), mm
\(F\# \ge\): f-numbers at or above this keep the ERM path error below \(\lambda/4\)

From \(z\tan^4\alpha/16 \le \lambda/4\) with \(\tan\alpha = 1/(2F\#)\). For θ ≠ 0 the fit point is not the apex of \(\tau'\) (which sits at \(u = z_0\tan\theta\)), and the error gains an odd, cubic part: at θ = 10°, \(z_0 = 30\) mm it is +69 µm at α = +18.4° and −24 µm at −18.4°. An odd phase error shifts the peak: 35 µm laterally at 9° in m24_stolt_planewave.m.

Demo 24.3

The ERM fit: true travel time against the fitted hyperbola

Exact formulaTravel times (24.8)–(24.9) and the fitted hyperbola from (24.11), evaluated as is.

true \(\tau'\) after compensationERM hyperbolaθ = 0 reference±λ/4

Steering θ 10°

Depth z₀ 30 mm

Centre frequency 5.0 MHz

Receive F#

Try: at θ = 0 the error is even and tiny inside F# 1.5 at 30 mm. Steer to 10°: the curve tilts and the error becomes lopsided (the cubic term). Go deep (80 mm) or high in frequency (10 MHz): the λ/4 band is crossed inside the aperture, and the readout's minimum F# rises. The exact map of Section 24.4 has zero error on this plot.

24.4The exact map: transmit plus receive wavevectors

No fit is needed. In the Born approximation and a homogeneous medium, write the incident plane wave as a delay \(e^{-j2\pi (f/c)(x\sin\theta + z\cos\theta)}\) and the receive leg as an angular spectrum of plane waves (the 2D Green's function is \(\int e^{\,j2\pi(k_x x + k_z\lvert z\rvert)}/k_z\, dk_x\) up to a constant). The x-integral over the object then picks out one object wavenumber per data sample:

\[ \begin{aligned} D(k_x, f) &\propto \frac{1}{k_z}\, O(K_x, K_z) \\ K_x &= k_x + (f/c)\sin\theta \\ K_z &= k_z + (f/c)\cos\theta \\ k_z &= \sqrt{(f/c)^{2} - k_x^{2}} \end{aligned} \]

(24.14)

\(O(K_x, K_z)\): 2D Fourier transform of the reflectivity
\((f/c)(\sin\theta, \cos\theta)\): transmit wavevector \(\mathbf{k}_{\mathrm{tx}}\)
\((k_x, k_z)\): receive wavevector \(\mathbf{k}_{\mathrm{rx}}\), length \(f/c\)
\(D(k_x, f)\): data spectrum, as in (24.1)
\(k_x,\ k_z,\ f,\ c,\ \theta\): as in (24.1) and (24.8)
\(K_x,\ K_z\): object wavenumbers, cycles/mm
\(1/k_z\): obliquity weight of the receive wave
\(\propto\): proportional to

Each sample measures the object at the transmit wavevector plus a receive wavevector of length \(f/c\). For fixed f these points lie on a circle of radius \(f/c\) centred on the transmit wavevector: the Ewald circle of diffraction tomography (Kak & Slaney, 1988; Figure 24.B). Lu's Fourier method (Cheng & Lu, 2006) maps plane-wave data the same way.

Schematic Figure 24.B · The Ewald circle of one plane wave: K = k_tx + k_rx; the receive fan selects an arc, the pulse band sweeps it into a crescent (to scale for θ = 15°, 5 MHz, F# 1.5; geometry from (24.14), not from data)

Inverting the map. After steering compensation the data column is \(K_x = k_x + (f/c)\sin\theta\), so again only f must be found. With \(K_x' = K_x - (f/c)\sin\theta\), squaring \(K_z - (f/c)\cos\theta = \sqrt{(f/c)^2 - K_x'^2}\) gives a closed form:

\[ \begin{aligned} f &= \frac{c\,(K_x^{2} + K_z^{2})}{2M} \\ M &= K_x\sin\theta + K_z\cos\theta \\ \left\lvert\frac{\partial f}{\partial K_z}\right\rvert &= \frac{c\,\lvert 2K_z M - (K_x^{2} + K_z^{2})\cos\theta\rvert}{2M^{2}} \end{aligned} \]

(24.15)

\(f\): data frequency to read for output sample \((K_x, K_z)\), MHz
\(M\): projection of \(\mathbf{K}\) on the transmit direction, cycles/mm
validity: \(K_z \ge (f/c)\cos\theta\) (upper half of the circle)
\(K_x,\ K_z,\ \theta,\ c\): as in (24.14)
\(\partial f/\partial K_z\): Jacobian of the change of variables, mm/µs
\(\lvert\cdot\rvert\): absolute value

θ = 0: \(f = c(K_x^2 + K_z^2)/(2K_z)\). For small \(K_x\), \(K_z \approx 2f/c - cK_x^2/(2f)\): exactly the second-order behaviour of the ERM, whose curve \(\sqrt{4(f/c)^2 - 2K_x^2}\) is an ellipse that osculates the Ewald circle at α = 0. They part at large receive angles, which is the ERM error of Section 24.3 seen in k-space.

The exact map costs the same as the ERM (one 1D interpolation per column) and removes the fit error: in m24_stolt_planewave.m its peak position error at 9° is below the 2 µm resolution of the measurement, against 35 µm for the ERM, with identical widths.

k-space support. One plane wave, one frequency and a receive half-angle \(\alpha_{\max}\) cover an arc of the Ewald circle. The pulse band sweeps the arc into a crescent; steering rotates the crescent's centre about the origin. For a plane-wave set:

\[ \begin{aligned} K_x &\ge (f/c)(\sin\theta_{\min} - \sin\alpha_{\max}) \\ K_x &\le (f/c)(\sin\theta_{\max} + \sin\alpha_{\max}) \\ K_z &\ge (f_{\mathrm{lo}}/c)(\cos\theta_{\max} + \cos\alpha_{\max}) \\ K_z &\le 2f_{\mathrm{hi}}/c \end{aligned} \]

(24.16)

\(\theta_{\min}, \theta_{\max}\): extreme steering angles of the set
\(\alpha_{\max}\): receive half-angle, \(\arctan(1/(2F\#))\)
\(f_{\mathrm{lo}}, f_{\mathrm{hi}}\): band edges, MHz
\(K_x,\ K_z\): object wavenumbers, cycles/mm, as in (24.14)
\(f,\ c\): as in (24.1); \(f/c\) is the wavenumber of the carrier, cycles/mm
\(\theta\): steering angle, as in (24.8)
\(\alpha\): receive arrival angle, as in (24.1)

The PSF is the inverse Fourier transform of the weighted coverage: Module 23's coarray, now with depth. Narrowband at 5 MHz, F# 1.5 (\(\sin\alpha = 0.316\)), a continuous span of ±9° gives a −6 dB width of 0.54 mm (receive-only: 0.59 mm, using \(\sin\alpha\) rather than the small-angle \(1.207\lambda F\# = 0.56\) mm).

Compounding in k-space. Images add, so spectra add. On a common f-k grid, sum the remapped spectra of all angles and apply one inverse 2D FFT per frame instead of one per angle. Where crescents overlap, spectra add coherently; gaps between them, from a coarse angle step, are the transmit replicas of Module 23.

Demo 24.4

k-space coverage of a plane-wave set and the PSF it implies

Approx. modelBorn, homogeneous medium: coverage from the exact map (24.14) with a Gaussian pulse spectrum and uniform receive \(k_x\) weight; PSF by a 256 × 256 inverse FFT of that coverage (no channel data, no element directivity).

Ewald circle at \(f_0\), \(\theta = \theta_{\max}\)ERM ellipse at \(f_0\), same θlateral PSF profile

Angles N 7

Span ±θ_max 9°

Receive F# 1.5

Fractional bandwidth 60 %

show ERM ellipse

Try: one angle gives a single crescent and the receive-only width. Open the span to ±20° with 21 angles: the crescents merge into a wide fan and the PSF narrows. Keep ±20° but drop to 3 angles: three separate crescents, and the PSF grows replica lobes. Raise F# to 3: each crescent thins. Steer far with the ERM ellipse shown: it hugs the circle near the middle and parts at the edges. Coverage weight: Gaussian pulse spectrum, uniform in receive \(k_x\); PSF from a 256 × 256 inverse FFT.

24.5f-k against DAS: same physics, different bookkeeping

Why they agree. Under the same model (Born, homogeneous, point elements) DAS and f-k both apply the adjoint of the forward operator: DAS sums the data along each pixel's delay curve; f-k evaluates the same integral in the Fourier domain. By stationary phase the two differ only by smooth weights in k-space: the Jacobian against DAS's element density, spreading and per-pixel normalization. Weights change amplitudes and sidelobe shapes, not where energy focuses.

Where they differ in practice:

Issue	DAS	f-k (Stolt, ERM or exact)
Receive aperture, F#	per-pixel element window z/F#	fan filter \(\lvert k_x\rvert \le (f/c)\sin\alpha_{\max}\) before mapping: a constant-F# dynamic aperture
Apodization	element weights	a window across the fan in \(k_x\)
Spatial aliasing	elements outside the window never enter	folded energy that lands inside the fan passes and is mapped to the wrong place
Interpolation	in time, on a smooth IQ envelope	in frequency; error grows with \(\Delta f\,\lvert t - t_c\rvert\)
Depth gain	unit-sum weights: flat for this lossless simulation	no per-pixel normalization: \(\propto\sqrt{z}\) relative to unit-sum DAS here (stationary-phase amplitude of the data spectrum; the Jacobian depends on angle only and cannot set a depth gain)
Medium	any delay law: layers, refraction, measured c maps	homogeneous (Stolt); horizontally layered by phase shift (Section 24.6)
Pixel grid	arbitrary	the FFT grid (x at pitch/U, z at the \(K_z\) step); other grids need a final interpolation
Edges	none	circular wrap in x and t unless zero-padded
Cost per angle	\(N_{\mathrm{pix}} \times N_{\mathrm{active}}\) interpolations	2D FFTs of the padded data, one interpolation per output sample, one \(\mathrm{IFFT}_2\)

Cost. DAS needs one interpolation per pixel per active element; f-k replaces the element factor by the logarithm of the FFT. For the 390 × 161-pixel grid of the MATLAB lab, DAS uses 3.4 million (pixel, element) pairs per angle; f-k uses FFTs of 4096 × 128 and 2 × 1721 × 512 points plus 0.39 million interpolations. Octave runs f-k 2.1–2.5× faster (0.29–0.39 s against 0.76–0.79 s per angle for two phantoms; the ratio varies from run to run with machine load). The ratio grows with pixel count and element count, and shrinks if DAS is restricted to a few lines or a small ROI. Demo 24.5 times both in your browser.

Demo 24.5

f-k against DAS on the same channel data: images, difference, timing

SimulationPlane-wave IQ channel data of points, speckle and a cyst synthesized in JS (128 elements, 0.30 mm, 5 MHz); f-k (exact or ERM map) and DAS run on the same data in the page and are timed there.

DASf-kdifference: f-k minus DAS envelope, ±6 dB

f-k map

Angles (3° step)

Interpolation · FFT length

Receive F# (DAS window, f-k fan)

Assumed c 1540 m/s

Dynamic range 50 dB

Try: the default images agree to the eye. The difference map is pink on average (f-k brighter with depth: the \(\sqrt{z}\) gain, since both are normalized on the 12 mm targets) and the speckle differs pixel by pixel by a few dB, because the k-space weights differ slightly; targets, sidelobes and the cyst agree. Pick "Nearest · 1×" or "Linear · 1× · not centred": haze fills the image at the depths of the targets and the speckle and the cyst contrast collapses from about −17 dB to about −2 dB, the interpolation smear at \(\varphi \approx 2\pi\). Pick 7 angles: sidelobes drop in both. Select the 9° steered set and switch to the ERM map: a small lateral offset appears in the readout. Set c to 1480 m/s: both methods defocus and shift alike, because both assume one speed. Above the targets f-k shows a faint haze where DAS is exactly zero (every echo from deeper scatterers arrives after the delays of a shallower pixel): the hard edge of the fan filter rings in both directions in z. A raised-cosine taper across the fan removes it (points only, 0°, in MATLAB: mean level above the targets −58 dB with a hard fan, −94 dB tapered). Timing is honest but not like for like: f-k computes the whole 512 × 512 FFT grid (262 k pixels), DAS only the 43 k displayed pixels with a precomputed delay table; compare the per-pixel cost. Point targets at 12 mm, anechoic cyst (radius 3 mm) in speckle at 24 mm; 128 elements, 0.30 mm pitch, 5 MHz.

The full comparison, simulated in MATLAB (m24_stolt_planewave.m: 7 angles over ±9°, F# 1.5, 0.15 × 0.077 mm pixels; widths and positions from spline-refined profiles):

Set	Method	Lateral −6 dB (mm)	Axial −6 dB (mm)	PSL (dB)	Peak offset (mm)	Cyst contrast (dB)	Beside block (dB)	Residual vs DAS (dB)
0°	DAS	0.535	0.234	−20.5	0	−16.6	−18.7	ref.
0°	f-k ERM	0.565	0.236	−20.8	0	−16.6	−18.1	−10.4
0°	f-k exact	0.555	0.236	−21.6	0	−16.8	−18.3	−11.6
7 angles	DAS	0.505	0.236	−35.1	0	−26.5	−28.3	ref.
7 angles	f-k ERM	0.505	0.238	−33.6	0	−26.9	−27.2	−16.9
7 angles	f-k exact	0.505	0.238	−34.2	0	−27.8	−27.6	−18.7
9° only	DAS / ERM / exact	0.555 / 0.555 / 0.550	0.236 / 0.238 / 0.238	−21.1 / −20.7 / −21.4	0 / 0.035 / 0	–	–	–

"Beside block" is the mean intensity in an empty strip 1.5 mm beside the speckle block relative to the block (clutter that leaked sideways). The residual is the complex difference after one complex gain per depth row, which removes the depth-gain difference: the f-k peak at 22 mm is 1.20× the one at 15 mm, against 0.95–0.99 for DAS, and \(\sqrt{22/15} = 1.21\). What remains is the k-space weighting. Interpolation (0°, exact map): nearest-bin on a 2048-point FFT gives contrast −8.2 dB and beside-block −7.2 dB; linear 2048 gives −11.3 and −10.6; linear 4096 gives −16.7 and −18.2, equal to DAS; centring or 8192 points change nothing further here.

Ultrasound link

Interpolation smear and aliasing leakage both look like clutter, and both fill cysts. If an f-k implementation shows lower contrast than DAS on the same data, check zero-padding and the fan filter before blaming the method: in the simulation above, contrast went from −11.3 dB to −16.7 dB (DAS −16.6 dB) by doubling the FFT length alone.

24.6Beyond one speed of sound: layers, sound-speed sweeps and the reflection matrix

Phase shift through layers. In a horizontally stratified medium \(c(z)\), a flat interface conserves \(k_x\) (Snell's law), so the angular spectrum can be stepped down slab by slab with the local speed. For a 0° plane wave the transmit leg is vertical and the step is

\[ \begin{aligned} \Psi(k_x, f; z + \Delta z) &= \Psi(k_x, f; z) \\ &\quad \times e^{\,j2\pi\left(k_z(z) + f/c(z)\right)\Delta z} \\ k_z(z) &= \sqrt{(f/c(z))^{2} - k_x^{2}} \\ \mathrm{image} &= \mathrm{IFFT}_{k_x} \textstyle\sum_f \Psi(k_x, f; z) \end{aligned} \]

(24.17)

\(\Psi(k_x, f; z)\): data spectrum continued down to depth z (receive leg plus vertical transmit leg)
\(c(z)\): local speed of the slab at depth z, mm/µs
\(\Delta z\): depth step (slab thickness), mm
\(k_x,\ f,\ z\): as in (24.1) and (24.3)
\(k_z(z)\): axial wavenumber of the receive wave in the slab, cycles/mm
\(f/c(z)\): axial wavenumber of the vertical transmit leg, cycles/mm
\(\sum_f\): sum over temporal frequencies
\(\mathrm{IFFT}_{k_x}\): inverse FFT over the lateral wavenumber \(k_x\), giving lateral position \(x\)
\(\mathrm{image}\): beamformed image line set at depth z
\(j\): imaginary unit

Gazdag (1978). Exact for horizontal layers (ignoring transmission losses and multiples); cost \(N_z N_f N_{kx}\) instead of one interpolation. Lateral variation needs split-step corrections (a phase screen per slab in x, then a homogeneous step; Stoffa et al., 1990) or a return to time-domain delays.

An 8 mm layer at 1450 m/s over tissue at 1540 m/s (m24_phase_shift_layered.m, ray-traced data): phase-shift migration puts every target within 1 µm of the truth with a 0.555 mm width at all depths. Stolt at 1540 m/s puts every target \((c_2/c_1 - 1)\,d_1 = 0.497\) mm too deep (measured 0.499 mm) and widens it by 15 %. Stolt with the time-average speed to 25 mm (1510 m/s) is exact only there: +0.20 mm at 15 mm, −0.20 mm at 35 mm.

Sound-speed sweeps. At θ = 0 the exact map is \(f = c\, g(K_x, K_z)\): changing c rescales the frequency at which each output sample reads the spectrum. A focus-versus-c sweep (Module 14's autofocus) therefore reuses the data FFTs and costs one interpolation and one IFFT per trial speed. The figure below shows how sharply coherence across the aperture peaks at the true c.

The reflection-matrix view. A full plane-wave (or element) acquisition is a reflection matrix R(input, output, f). f-k imaging propagates both sides to each depth and keeps the confocal part, input and output focused at the same point. The distortion-matrix approach (Lambert et al., 2020) keeps the off-confocal terms to estimate and correct aberrations. REFoCUS (Bottenus, 2018) recovers the element-wise synthetic-aperture data set from focused transmits by undoing the known transmit delays frequency by frequency, after which any of these Fourier or matrix methods applies.

Name clash

"Fourier-domain beamforming" also names a different idea: Chernyakova & Eldar (2014) compute each scanline's DAS output from a small set of Fourier coefficients of the channel data, to sample below Nyquist. That is per-line DAS in the frequency domain, not migration; it keeps DAS's geometry and cost structure.

24.7Visual summary

Fourteen pictures that carry the module. Each card is a question to answer in five seconds, the answer with its number, and a figure computed from the formula; the badge on each card says whether it is a closed form or a simplified model.

24.8How experts think about it

An expert treats an imaging algorithm as a map from data wavenumbers \((k_x, f)\) to object wavenumbers \((K_x, K_z)\). The physics fixes the map (transmit plus receive wavevector); an algorithm is good to the extent it reproduces the map's phase, and the amplitudes are a choice of weighting. Every failure then has an address: a wrong phase (ERM fit, wrong c, layers), a wrong sample (interpolation, aliasing, wrap), or a weight (Jacobian, fan, apodization).

Every point of the channel-data spectrum is a plane wave: \(k_x = (f/c)\sin\alpha\).

cone \(\lvert k_x\rvert \le f/c\)5 MHz: 3.25 /mmα = 30° → 1.62 /mm

A plane wave crossing the array at angle α has a linear delay \(x\sin\alpha/c\), a phase ramp \(2\pi f x\sin\alpha/c\), so its lateral wavenumber is \((f/c)\sin\alpha\). Straight lines through the origin of the \((k_x, f)\) plane are arrival directions; the cone edge is 90°. For pulse-echo plane-wave data the receive leg is one-way, so the data fill the cone; reading an angle off a spectral feature is one division.

Exact formula

Homogeneous medium under the array; a lens or matching layer refracts and rescales angles; zero-offset data fill \(2f/c\).

Outside the cone there is nothing to image, and nothing to back-propagate: \(1.1\,f/c\) decays 81 dB per millimetre.

\(e^{-2\pi z\sqrt{k_x^2 - (f/c)^2}}\)\(1.01\,f/c\): 25 dB/mm\(1.5\,f/c\): 198 dB/mm

Beyond \(\lvert k_x\rvert = f/c\) the wave equation makes \(k_z\) imaginary. Scatterers are millimetres away, so their evanescent parts arrive more than 100 dB down. Run the propagator backwards and the decay becomes a gain of the same size, applied to noise and leakage, which is why every Fourier imager masks the region (Stolt never reads it: \(f \ge v\lvert k_x\rvert\)). Energy you see outside the cone is aperture-truncation leakage: −22 dB for a 5 mm target with a rectangular aperture.

Exact formula

Near-field scatterers within a fraction of λ (lens, skin line); leaky or surface waves in the probe; truncation leakage is not evanescent energy.

Plane-wave receive needs λ/2 pitch, zero-offset data need λ/4; at λ pitch every arrival steeper than 31° folds at 5 MHz.

PW: \(p \le \lambda/(2\sin\alpha)\)ZO: \(p \le \lambda/(4\sin\alpha)\)0.30 mm: \(\alpha_{\mathrm{alias}}\) 30.9° (5 MHz), 20.0° (7.5 MHz)

The data cone reaches \(k_x = (f/c)\sin\alpha\) (one-way) or \(2(f/c)\sin\alpha\) (two-way); Nyquist is \(1/(2p)\). Setting them equal gives the limits. For a 0° plane wave and 0.30 mm pitch, 82 % of the channel energy of a target at 5 mm lies beyond Nyquist, 10 % at 30 mm. With a receive F# the angle is bounded (\(\sin\alpha = 0.316\) at F# 1.5), so the fan itself stays inside Nyquist for \(p \le 1.58\lambda\) (plane waves) and \(0.79\lambda\) (zero-offset). That does not make the fan alias-free: steep arrivals outside the fan fold by \(1/p\) into it unless \(p \le \lambda/(1 + \sin\alpha) = 0.76\lambda\) (plane waves) or \(\lambda/(2(1 + \sin\alpha)) = 0.38\lambda\) (zero-offset), or element directivity suppresses them.

Exact formula

Element directivity attenuates steep arrivals but does not remove them; aliasing that has folded cannot be filtered afterwards.

Pulse-echo is one-way propagation at \(c/2\): migrate zero-offset data with half the speed of sound.

\(\tau = 2r/c = r/(c/2)\)\(v = c\) → depth ×210 mm → 19.94 mm, −21 dB

The exploding-reflector model (24.4): zero-offset delays equal one-way delays from sources at the reflectors that fire at t = 0 in a medium of speed \(c/2\). Migrating with \(v = c\) halves every \(k_z\), so depths double and the hyperbola curvatures no longer match (m24_fk_migration_zero_offset.m: 21–24 dB peak loss). For plane waves the equivalent speed is Garcia's \(\hat{v}\), \(c/\sqrt{2}\) at broadside.

Exact formula

Exact for zero-offset (monostatic) data only; any separation of transmit and receive breaks the one-way equivalence.

Stolt migration is time-to-depth conversion plus a dip correction: \(k_z = \sqrt{(2f/c)^2 - k_x^2}\).

\(k_x = 0\): \(k_z = 2f/c\)steep \(k_x\): largest shift

Each sample moves straight down its \(k_x\) column from \(2f/c\) to the circle of radius \(2f/c\) (Figure 24.A, (24.6)). Flat reflectors (\(k_x = 0\)) are only relabelled, which is time-to-depth conversion. A sample at \(k_x = 0.8 \times 2f/c\) moves to \(0.6 \times 2f/c\), the largest correction: steeply dipping energy, the hyperbola flanks, collapses most. Total cost: two FFTs and one interpolation per column.

Exact formula

Homogeneous medium, one speed; the output grid is the FFT grid.

The Jacobian is an obliquity weight, the cosine of the propagation angle; it changes amplitudes, not focus.

\(v k_z/\lvert\mathbf{k}\rvert = v\cos\varphi\)30°: −1.2 dB60°: −6.0 dB

\(df = \lvert\partial f/\partial k_z\rvert\, dk_z\) converts the frequency integral into a \(k_z\) integral; for zero-offset Stolt it is \(v k_z/\sqrt{k_x^2 + k_z^2} = v\cos\varphi\). It is real and smooth, so it reshapes the PSF's k-space weighting slightly but never moves energy. Being a weight on \((K_x, K_z)\) alone, it acts as a shift-invariant filter and cannot create a depth-dependent gain. Demo 24.2 lets you switch it off: the image barely changes. Phase-shift migration needs no weight at all: it integrates over f directly (m24_kspace_evanescent.m: Stolt and phase-shift peaks within 2.1 %, phases equal).

Exact formula

For true-amplitude work (attenuation, backscatter) the weights matter and must be derived for the actual source and receiver model.

Interpolation along f is Stolt's only error in a homogeneous medium: keep \(2\pi\,\Delta f\,\lvert t - t_c\rvert\) below about one.

smear \(0.037\,\varphi^2\)φ = 0.92 → −30 dBφ = 3 → −10 dB

Linear interpolation of \(e^{-j2\pi f t}\) at fractional bin w errs by \(w(1 - w)\varphi^2/2\); the column-to-column fluctuation (rms \(0.037\,\varphi^2\)) smears each echo across x at its depth (24.7). Zero-padding shrinks Δf, centring the time origin shrinks \(\lvert t - t_c\rvert\) by half. The figure plots the smear at the end of the record against the padding factor. Simulated: 2048 points, uncentred, 62.5 µs record → beside-block clutter −10.6 dB (theory at φ = 3: −10.0 dB); 4096 points → −18.2 dB, equal to DAS.

Approx. modelSmear law \(0.037\,\varphi^2\): coefficient from a Monte Carlo over fractional bin positions, linear interpolation, one echo.

Higher-order or windowed-sinc interpolators cut the error further at more cost; nearest-bin is never acceptable.

Steering compensation is a spectral shear: multiply by \(e^{\,j2\pi f x\sin\theta/c}\) and the cone tilts by \(\sin\theta\).

\(k_x \to k_x + f\sin\theta/c\)15°, 7 MHz: +1.18 /mm

Removing the transmit delay \(x\sin\theta/c\) from each channel makes the plane-wave delay curve depend only on \(x - x_0\), the shift invariance a migration needs (24.9). In k-space it adds \(f\sin\theta/c\) to every \(k_x\): a shear growing with frequency. The sheared spectrum can cross the lateral Nyquist even when the raw data do not, so upsample in x (or apply the shear inside the mapping) before transforming back.

Exact formula

Assumes the wavefront passes the array centre at t = 0; any other time origin is a further phase \(e^{\,j2\pi f t_0}\).

Garcia's ERM speed is \(c/\sqrt{1 + \cos\theta + \sin^2\theta}\): \(c/\sqrt{2}\) at broadside.

θ = 0: \(\hat{v} = 0.7071c\), \(\beta = \sqrt{2}\)10°: 0.7045c, 1.388, γ 0.17120°: 0.6973c, 1.314, 0.323

Matching the compensated travel time \(z_0\cos\theta - u\sin\theta + \sqrt{u^2 + z_0^2}\) by a hyperbola \(\sqrt{(u - \hat{u})^2 + \hat{z}^2}/\hat{v}\) in value, slope and curvature at u = 0 gives three equations and three unknowns: \(\hat{v}\), \(\hat{z} = \beta z_0\) and \(\hat{u} = \gamma z_0\) (24.11). At θ = 0, apex \(2z_0/c\) and curvature \(1/(z_0 c)\) force \(\hat{z} = \sqrt{2}\,z_0\) and \(\hat{v} = c/\sqrt{2}\). The image is mapped back by \(\hat{k}_z = (K_z - \gamma K_x)/\beta\).

Exact formula

A local fit: exact to second order at the element above the scatterer only.

The ERM fit errs by \(z\tan^4\alpha/16\); f-k with the ERM is safe for \(F\# \gtrsim (z/64\lambda)^{1/4}\), about 1.1 at 30 mm.

30 mm, 18.4°: 21 µm26.6°: 99 µm\(F\#_{\min}\) 1.00 / 1.11 / 1.32 at 20 / 30 / 60 mm

Expand true and fitted travel times at θ = 0 with \(t = \tan\alpha\): \(z(2 + t^2/2 - t^4/8)\) against \(z(2 + t^2/2 - t^4/16)\). The difference \(z t^4/16\) must stay below λ/4; with \(\tan\alpha = 1/(2F\#)\) that gives \(F\# \ge (z/(64\lambda))^{1/4}\) (24.13). The bound rises slowly with depth and frequency. Steered waves add an odd error that shifts the peak (35 µm at 9°). The exact map has none of this at equal cost.

Exact formula

Leading-order term: overestimates the exact error at large angles (117 vs 99 µm at 26.6°); λ/4 is a convention, not a hard limit.

Each data sample measures the object at transmit plus receive wavevector: an Ewald circle of radius \(f/c\).

\(f = c\lvert\mathbf{K}\rvert^2/(2\,\mathbf{K}\cdot\hat{\mathbf{n}}_\theta)\)θ = 0: \(K_z = f/c + \sqrt{(f/c)^2 - k_x^2}\)

Born plus the angular-spectrum Green's function give \(D \propto O(\mathbf{k}_{\mathrm{tx}} + \mathbf{k}_{\mathrm{rx}})/k_z\) (24.14). Inverting for f is closed form (24.15), so the exact map costs one interpolation per column, like Stolt. The ERM curve (dashed) is an ellipse osculating the circle at α = 0; at θ = 0 both expand as \(K_z \approx 2f/c - cK_x^2/(2f)\).

Exact formula

Born approximation, homogeneous medium, infinite plane wave; finite apertures add edge waves the map does not model.

In f-k, compound before the inverse FFT: spectra add, so N angles cost one \(\mathrm{IFFT}_2\).

\(\sum_\theta \mathrm{IFFT}(R_\theta) = \mathrm{IFFT}(\sum_\theta R_\theta)\)crescent centres at \((f/c)\,\hat{\mathbf{n}}_\theta\)

Coherent compounding sums complex images; the inverse FFT is linear, so summing remapped spectra on a common grid and transforming once is identical. Each angle contributes a crescent centred on its transmit wavevector; the span sets the lateral extent of the union (resolution), overlaps add coherently, and gaps from a coarse step are the replicas of Module 23. The figure draws the −6 dB band crescents for −12° to +12° in 6° steps at F# 1.5.

Exact formula

Motion between angles makes the spectra add with phase errors (Module 23's motion rule applies unchanged).

A k-space fan filter is a constant-F# dynamic aperture, and it cannot undo aliasing that has already folded into the fan.

\(\lvert k_x\rvert \le (f/c)\sin(\arctan(1/2F\#))\)F# 1.5: \(0.316\,f/c\)

By stationary phase, receive angle α maps to \(k_x = (f/c)\sin\alpha\) at every depth, so masking the fan applies the same angular aperture everywhere, which is what a constant-F# DAS window does until the aperture saturates. But a steep arrival beyond Nyquist folds by \(1/p\) and can land inside the fan (red line: the cone folded for 0.30 mm pitch). DAS never uses those elements; f-k has already mixed them in.

Exact formula

Near the array edges DAS's aperture saturates and f-k's does not; a hard fan edge rings in z (taper it, the k-space version of apodization); apodization inside the fan is a \(k_x\) window, not an element window.

f-k and DAS differ by smooth k-space weights: expect equal widths, a \(\sqrt{z}\) gain difference and a −17 to −19 dB residual.

width 0.505 = 0.505 mm (7 angles)A(22)/A(15): 1.20 vs \(\sqrt{22/15} = 1.21\)

Both are adjoints of the same Born operator, so they focus identically; stationary phase leaves only amplitude factors (Jacobian, element density, spreading, DAS's per-pixel normalization). Simulated with 7 angles: lateral width equal, PSL within 1.5 dB, cyst contrast −26.5 (DAS) against −26.9 and −27.8 dB, and a residual of −17 to −19 dB once one gain per depth row removes the \(\sqrt{z}\) trend shown here.

Approx. modelLine: \(\sqrt{z/15}\) depth-gain law of unnormalized f-k (stationary-phase amplitude of a constant-amplitude echo's spectrum) against flat unit-sum DAS; dots: the two simulated values at 22 mm (m24_stolt_planewave.m).

Holds only if both use the same aperture (fan versus window), the same c, and f-k is padded and fan-filtered.

DAS costs pixels × elements, f-k costs N log N: the speed-up is of order \(N_{\mathrm{active}}/\log_2 N\).

Octave: 2.1–2.5× (63 k pixels, 54 active)flop model: 2.2×DAS 3.4 M pairs / angle

DAS does an interpolation and a complex multiply per (pixel, active element), about 14 flops. f-k does the equivalent of two 2D FFTs on a 2× padded grid (\(5N\log_2 N\) flops each) plus one interpolation per output sample. With 128 elements at F# 1.5 that predicts 2.2× for the lab's 390 × 161 grid; Octave measured 2.1–2.5×. The gain is modest for a 1D array and grows with the element count: the dashed curve is a 500-element active aperture, typical of a 2D matrix array, where Fourier methods pay most.

Approx. modelFlop counts: DAS 14 flops per (pixel, element); f-k \(20N\log_2(2N) + 10N\); memory traffic ignored.

Small ROIs or few scanlines favour DAS; f-k's grid is fixed by the FFT; memory traffic, not flops, often dominates, and a GPU DAS narrows the gap.

Stolt assumes one speed: under a layer every target shifts by \((c_2/c_1 - 1)\,d_1\), and phase-shift migration removes it.

8 mm at 1450 m/s: +0.497 mmmeasured +0.499 mm\(\bar{c}\) fixes one depth only: ±0.20 mm

The two-way time to a target below the layer is preserved; a single speed \(c_2\) converts the extra time \(2d_1(1/c_1 - 1/c_2)\) into depth \(c_2 d_1(1/c_1 - 1/c_2)\). An average speed \(\bar{c} = z_{\mathrm{ref}}/t_{\mathrm{ref}}\) zeroes the error at \(z_{\mathrm{ref}}\) only (lines: 1540 m/s flat, 1510 m/s crossing zero at 25 mm). Phase shift with the local speed conserves \(k_x\) across the interface and focuses everywhere (24.17) (width 0.555 mm at every depth, against 0.635 mm for Stolt at 1540 m/s).

Approx. modelVertical rays, horizontal layer, depth = assumed speed × one-way time; refraction of oblique rays ignored.

Horizontal layers only; tilted or rough interfaces need split-step or time-domain delays; transmission losses and multiples ignored.

A sound-speed sweep in f-k reuses the FFTs: the map \(f = c\,g(\mathbf{K})\) only rescales the frequency axis.

θ = 0: \(f = c(K_x^2 + K_z^2)/(2K_z)\)per trial c: 1 interpolation + 1 IFFT

At broadside the exact (and ERM) maps are linear in c: every output sample reads the same data spectrum at a frequency proportional to c. A focus metric versus c (coherence, sharpness; Module 14) needs no new forward FFTs, which makes dense sweeps cheap. The figure draws \(f(K_z)\) at \(K_x = 2\) /mm for 1400, 1540 and 1700 m/s. Steered angles add c to the steering compensation and the fan, still without new data FFTs in time.

Exact formula

A global c only; local speed estimation needs layered or local methods (Module 14, Section 24.6).

24.9Cheat sheet

Concept	Formula	Ultrasound use
Plane-wave component	\(k_x = \frac{f}{c}\sin\alpha,\ \ k_z = \sqrt{(f/c)^2 - k_x^2}\)	read arrival angles off the data spectrum
Evanescent decay	\(e^{-2\pi z\sqrt{k_x^2 - (f/c)^2}}\); 81 dB/mm at \(1.1\,f/c\), 5 MHz	mask outside the cone; never back-propagate
Phase-shift propagator	\(\times\, e^{\,j2\pi k_z\Delta z}\)	angular spectrum method, layered migration
Spatial Nyquist	\(\lvert k_x\rvert \le \frac{1}{2p}\); PW \(p \le \frac{\lambda}{2\sin\alpha}\), ZO \(p \le \frac{\lambda}{4\sin\alpha}\)	pitch choice, alias angle \(\arcsin(c/(2pf))\)
Exploding reflector	\(\tau = \frac{2r}{c} = \frac{r}{c/2}\)	zero-offset migration speed \(v = c/2\)
Stolt (zero-offset)	\(f = \frac{c}{2}\sqrt{K_x^2 + K_z^2},\ \ J = \frac{c}{2}\frac{K_z}{\lvert\mathbf{K}\rvert}\)	SAFT imaging in N log N
Steering compensation	\(\times\, e^{\,j2\pi f x\sin\theta/c} \Leftrightarrow k_x \to k_x + \frac{f\sin\theta}{c}\)	first step of every plane-wave f-k
Garcia ERM	\(\hat{v} = \frac{c}{\sqrt{1 + \cos\theta + \sin^2\theta}},\ \beta = \frac{(1 + \cos\theta)^{3/2}}{1 + \cos\theta + \sin^2\theta},\ \gamma = \frac{\sin\theta}{2 - \cos\theta}\)	plane-wave Stolt; map \(f = \hat{v}\sqrt{K_x^2 + ((K_z - \gamma K_x)/\beta)^2}\)
ERM error	\(z\tan^4\alpha/16\) (θ = 0); \(F\# \ge (z/(64\lambda))^{1/4}\)	usable aperture of the ERM
Exact plane-wave map	\(f = \frac{c\,(K_x^2 + K_z^2)}{2\,(K_x\sin\theta + K_z\cos\theta)}\)	fit-free f-k, same cost
Ewald circle	centre \(\frac{f}{c}(\sin\theta, \cos\theta)\), radius \(f/c\)	k-space support, compounding, PSF prediction
Interpolation smear	\(0.037\,\varphi^2,\ \ \varphi = 2\pi\,\Delta f\,\lvert t - t_c\rvert\)	zero-padding and time-centring choice
Fan filter	\(\lvert k_x\rvert \le \frac{f}{c}\sin\!\big(\arctan\frac{1}{2F\#}\big)\)	F-number in f-k
Layer depth shift	\(\Delta z = \left(\frac{c_2}{c_1} - 1\right) d_1\)	fat layer error; phase-shift fix
Cost	DAS \(N_{\mathrm{pix}} N_{\mathrm{active}}\); f-k \(\approx N\log N + N_z N_x\)	real-time and GPU budgeting

24.10Interview and qualifying-exam questions

I1Why does Stolt migration of pulse-echo data use \(c/2\)? · screen

Zero-offset delays are \(2r/c\), identical to one-way delays \(r/v\) with \(v = c/2\) from sources placed at the reflectors and fired at t = 0 (exploding reflector, (24.4)). The data are then a one-way field, and migration back-propagates it to t = 0 with speed \(c/2\). Using c doubles depths (10 → 19.94 mm in the lab) and defocuses by about 21 dB. Follow-up: "And for plane waves?" There is no exact equivalent speed; Garcia's fit gives \(\hat{v} = c/\sqrt{1 + \cos\theta + \sin^2\theta}\), \(c/\sqrt{2}\) at 0°.

I2What is the evanescent region of a \((k_x, f)\) spectrum, and why is it empty in channel data? · screen

\(\lvert k_x\rvert > f/c\), where \(k_z\) is imaginary and the field decays as \(e^{-2\pi z\sqrt{k_x^2 - (f/c)^2}}\): 81 dB per mm at \(1.1\,f/c\) at 5 MHz. Scatterers are millimetres away, so nothing evanescent reaches the array; what shows up there is truncation leakage (−22 dB for a 5 mm target, rectangular aperture). Follow-up: "What if you back-propagate it?" The decay becomes a gain on noise; mask it.

I3Derive Stolt's mapping and its Jacobian. What happens if you drop the Jacobian? · onsite

\(r(x, z) = \iint D(k_x, f)\, e^{\,j2\pi(k_x x + k_z z)}\, df\, dk_x\) with \(k_z = \sqrt{(f/v)^2 - k_x^2}\). Substitute \(f = v\sqrt{k_x^2 + k_z^2}\), \(df = v k_z/\sqrt{k_x^2 + k_z^2}\, dk_z\); the image is an inverse 2D FFT of D resampled at that f times the Jacobian (24.6). Dropping it changes the k-space weighting (a \(\cos\varphi\) obliquity), hence slightly the PSF shape and amplitude, not the focusing and not the \(\sqrt{z}\) depth gain of unnormalized f-k. Follow-up: "Where does Stolt need interpolation and how fine?" Along f; keep \(2\pi\,\Delta f\,\lvert t - t_c\rvert \lesssim 1\).

I4Why can't zero-offset Stolt be applied to plane-wave data directly? Derive Garcia's ERM velocity. · onsite

The plane-wave delay \([x_0\sin\theta + z_0\cos\theta + \sqrt{(x - x_0)^2 + z_0^2}]/c\) is not a hyperbola and its transmit part depends on \(x_0\). Steering compensation (shift by \(-x\sin\theta/c\)) makes it depend on \(x - x_0\) only; then fit \(\sqrt{(u - \hat{u})^2 + \hat{z}^2}/\hat{v}\) in value, slope and curvature at u = 0. That gives \(R^2 = z_0^2(1 + \cos\theta)^2/(1 + \cos\theta + \sin^2\theta)\), \(\hat{v} = c/\sqrt{1 + \cos\theta + \sin^2\theta}\), \(\beta = (1 + \cos\theta)^{3/2}/(1 + \cos\theta + \sin^2\theta)\), \(\gamma = \sin\theta/(2 - \cos\theta)\) (24.10)–(24.11). Follow-up: "How big is the error?" \(z\tan^4\alpha/16\) at 0°: F# ≳ 1.1 at 30 mm, 5 MHz.

I5How do you implement an F-number and apodization in f-k? · onsite

Receive angle α corresponds to data wavenumber \(k_x = (f/c)\sin\alpha\), so the F# is a fan filter \(\lvert k_x\rvert \le (f/c)\sin(\arctan(1/(2F\#)))\) applied before steering compensation and mapping, and apodization is a window across the fan. It equals a constant-F# dynamic aperture. Follow-up: "Is it identical to DAS?" Not at the array edges (DAS's aperture saturates) and not for aliased data, which can fold into the fan.

I6Your f-k images show a lateral haze at the depth of every bright reflector, absent in DAS. Diagnose. · onsite

Interpolation along f is too coarse: the error \(w(1 - w)\varphi^2/2\) with \(\varphi = 2\pi\,\Delta f\, t\) varies from column to column and spreads each echo across x at its depth (rms \(0.037\,\varphi^2\); −10 dB at φ = 3). Fix: zero-pad the time FFT (2–4× the record) and centre the time origin. In the lab, doubling the FFT length took beside-block clutter from −10.6 to −18.2 dB (DAS −18.7). Second suspect: lateral wrap or aliasing; check padding in x and the fan. Follow-up: "Why only at the depth of the echo?" The error multiplies that echo's own spectrum.

I7Compare the cost of f-k and DAS. When is f-k not worth it? · onsite

DAS: \(N_{\mathrm{pix}} \times N_{\mathrm{active}}\) interpolations per angle (3.4 M in the lab). f-k: FFTs of the padded data (N log N), one interpolation per output sample, one IFFT; compounding adds spectra and needs one IFFT per frame. Measured 2.1–2.5× in Octave; the asymptotic gain is about \(N_{\mathrm{active}}/\log_2 N\). f-k loses for small ROIs or a few lines, arbitrary pixel grids, heterogeneous media, and when a GPU DAS is memory-bound anyway. Follow-up: "And for 3D with a matrix array?" The element factor grows to thousands, which is where Fourier methods pay most.

I8The exact plane-wave map is closed form. Why would anyone use the ERM fit? · research

Both are one 1D interpolation per column after the same steering compensation, so cost is not the reason. The ERM keeps the Stolt square-root structure, shares code with zero-offset migration, and its error is negligible at moderate F# and steering. The exact map (24.15) removes the fit error (35 µm lateral offset at 9° in the lab, none for the exact map). A defensible choice today: the exact map, unless the code base is built around ERM. Follow-up: "What does the exact map still assume?" Born, homogeneous c, an infinite plane wave and point elements.

I9How would you extend f-k imaging to a subcutaneous fat layer, and then to lateral heterogeneity? · research

Horizontal layers: phase-shift migration, stepping Ψ by \(e^{\,j2\pi(k_z(z) + f/c(z))\Delta z}\) (24.17); \(k_x\) is conserved by Snell's law, so it is exact. An 8 mm layer at 1450 m/s shifts targets by 0.497 mm under single-speed Stolt; phase shift removes it. Lateral variation: split-step Fourier (a phase screen \(e^{\,j2\pi f\Delta s(x)\Delta z}\) per slab, then a homogeneous step) or phase shift plus interpolation between reference speeds; beyond that, time-domain delays from an eikonal solver. Follow-up: "Where does the speed map come from?" From a focus or coherence sweep (cheap in f-k), or tomography (Module 13).

I10How can Fourier-domain imaging help sound-speed estimation and aberration correction research? · research

At broadside the map is \(f = c\,g(\mathbf{K})\), so each trial c reuses the data FFTs: a dense global-c sweep costs one interpolation and one IFFT per value. For local aberrations, keep the full reflection matrix (all transmit angles and receive elements) instead of the confocal sum: propagating both sides to a depth in the angular-spectrum domain gives focused-basis matrices whose off-confocal structure reveals the aberration (distortion-matrix approach). Follow-up: "And for focused-beam data?" REFoCUS recovers the synthetic-aperture data set first.

I11A 128-element probe at 0.30 mm pitch runs at 7.5 MHz. Will f-k images alias, and what do you do? · screen

Nyquist is \(1/(2 \times 0.30) = 1.67\) /mm; the receive cone at 7.5 MHz reaches 4.87 /mm, so arrivals steeper than \(\arcsin(c/(2pf)) = 20.0°\) fold. With F# 1.5 the fan needs \(0.316 \times 4.87 = 1.54\) /mm, still inside Nyquist, but steep energy folded from beyond can land in the fan. Use the fan filter, let element directivity help, and limit the angles that matter, or accept DAS where the aperture window keeps steep elements out. Follow-up: "Zero-offset at the same pitch?" Cone \(2f/c\), hopeless without a small aperture.

24.11MATLAB lab

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

Download the scripts: m24_stolt_planewave.m · m24_fk_migration_zero_offset.m · m24_kspace_evanescent.m · m24_phase_shift_layered.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m24_stolt_planewave.mPlane-wave Stolt (Garcia ERM and exact map) against DAS: timing, widths, PSL, positions, cyst, interpolation · change thetaDeg, Fnum, ntFFT, tCenter, useFan · runtime ≈ 15 s

%% m24_stolt_planewave.m  --  Stolt f-k migration for plane waves versus DAS
% Purpose : simulate complex baseband (IQ) channel data of a 128-element linear array for
%           steered plane waves (point targets, and an anechoic cyst in speckle), then form
%           the image three ways on the same pixel grid:
%             DAS   : delay and sum, dynamic receive aperture F# 1.5 (rectangular)
%             ERM   : Stolt f-k migration with the exploding-reflector fit of Garcia et al.
%                     (2013): v = c/sqrt(1+cos(th)+sin(th)^2), depth factor
%                     beta = (1+cos th)^1.5/(1+cos th+sin^2 th), lateral shear sin th/(2-cos th)
%             exact : the same pipeline with the exact transmit + receive wavevector map
%                     f = c(Kx^2+Kz^2) / (2(Kx sin th + Kz cos th))
%           Pipeline (both maps): FFT in t, FFT in x, receive fan filter equivalent to F# 1.5,
%           2x lateral upsampling, steering compensation exp(j2 pi f x sin(th)/c), ONE linear
%           interpolation along f per Kx column, Jacobian |df/dKz|, IFFT2.
%           Printed: time per angle, -6 dB widths, peak sidelobe, peak position error at
%           0 and 9 deg, cyst contrast, clutter beside the speckle block, residual against
%           DAS, and how interpolation, zero-padding and time-centering change the result.
% Change  : thetaDeg (compounding set), Fnum, ntFFT, tCenter, useFan, pitch.
% Expect  : widths equal to DAS within ~2 % (lateral 0.55 mm single angle, 0.50 mm for
%           7 angles; axial 0.24 mm), PSL within ~1.5 dB, cyst contrast within ~1 dB.
%           ERM at 9 deg: peak displaced laterally by tens of um; exact map: < 5 um.
%           Linear interpolation on a 2048-point FFT (no centering) smears every strong
%           echo laterally (clutter beside the block -11 dB instead of -19 dB); 4096 points
%           with time-centering removes it. f-k is a few times faster than DAS in Octave.
% Model   : linear acoustics, first-order Born, infinite plane-wave transmit, point-like
%           elements (no directivity), no attenuation or spreading, Gaussian two-way pulse
%           (60 % band). Runtime ~20 s in Octave.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);
tStart = tic;

%% Parameters
c     = 1.540;                 % speed of sound (mm/us)
f0    = 5;                     % center frequency (MHz)
lam   = c/f0;                  % wavelength (mm) = 0.308
FBW   = 0.60;                  % -6 dB fractional bandwidth (amplitude spectrum)
fs    = 20;                    % IQ sampling rate (MHz)
Nel   = 128; pitch = 0.30;     % elements, pitch (mm)
xe    = ((0:Nel-1) - (Nel-1)/2) * pitch;   % element positions (mm)
Fnum  = 1.5;                   % receive F-number (DAS aperture z/Fnum, f-k fan filter)
alphaMax = atan(1/(2*Fnum));   % receive half-angle allowed by the F-number (rad)
useFan = true;                 % receive fan filter in f-k (the k-space F-number)
thetaDeg = -9:3:9;             % plane-wave angles (deg); must contain 0 and 9
theta = thetaDeg*pi/180;  Nang = numel(theta);  iMid = find(thetaDeg == 0);
density = 40;                  % speckle scatterers per mm^2
nt    = 1250;                  % samples per channel (62.5 us), t = 0 at transmit
ntFFT = 4096;                  % time FFT length for f-k (3.3x zero-padding)
tCenter = nt/(2*fs);           % time-centering reference (us) before f interpolation

sig_f = FBW*f0/(2*sqrt(2*log(2)));          % Gaussian spectrum std (MHz) = 1.27
sig_t = 1/(2*pi*sig_f);                     % envelope std (us) = 0.125
kt    = (-ceil(4*sig_t*fs):ceil(4*sig_t*fs)).';
envK  = exp(-(kt/fs).^2/(2*sig_t^2));       % sampled envelope

%% Phantoms: (A) point targets, (B) speckle with an anechoic cyst
ptsX = [-10 0 10 0].';  ptsZ = [15 15 15 22].';  ptsA = ones(4,1);
R = 3; cx = 0; cz = 33;  xr = [-8 8]; zr = [27 39];
Ns  = round(density*diff(xr)*diff(zr));
spX = xr(1) + diff(xr)*rand(Ns,1);  spZ = zr(1) + diff(zr)*rand(Ns,1);
keep = (spX-cx).^2 + (spZ-cz).^2 > R^2;  spX = spX(keep);  spZ = spZ(keep);
spA = (randn(numel(spX),1) + 1i*randn(numel(spX),1))/sqrt(2);
fprintf('speckle scatterers %d, cyst radius %.0f mm, %d angles (%s deg)\n', numel(spX), R, Nang, mat2str(thetaDeg));

%% Image grids (fixed): lateral 2x upsampled, depth from a fixed Kz grid
nxp  = 256;  U = 2;  nxf = U*nxp;  dxf = pitch/U;
xf   = xe(1) + (0:nxf-1)*dxf;               % image x grid (mm), 0.15 mm
kxd  = [0:nxp/2-1, -nxp/2:-1]/(nxp*pitch);  % data wavenumbers (cycles/mm)
Kx   = [0:nxf/2-1, -nxf/2:-1]/(nxf*dxf);    % image wavenumbers (cycles/mm)
NZ   = 1024;  dKz = 2*(fs/2048)/c;          % Kz step (cycles/mm): depth extent 78.8 mm
Kz   = (0:NZ-1).'*dKz;  zf = (0:NZ-1).'/(NZ*dKz);   % image z grid (mm), 0.077 mm
[KX, KZ] = meshgrid(Kx, Kz);
ix = find(abs(xf) <= 12);  iz = find(zf >= 10 & zf <= 40);
xg = xf(ix);  zg = zf(iz);  [XG, ZG] = meshgrid(xg, zg);
fprintf('pixel grid %d x %d (dx %.3f mm, dz %.4f mm)\n', numel(zg), numel(xg), dxf, zg(2)-zg(1));

%% DAS geometry (angle independent)
Px = XG(:);  Pz = ZG(:);  nPix = numel(Px);
Rrx = sqrt((Px - xe).^2 + Pz.^2)/c;          % receive time (us), nPix x Nel
act = abs(xe - Px) <= Pz/(2*Fnum);           % dynamic aperture
[pp, ee] = find(act);  rxT = Rrx(act);  clear Rrx
Wn  = 1 ./ accumarray(pp, 1, [nPix 1]);      % unit-sum rectangular weights
Wc  = Wn(pp) .* exp(1i*2*pi*f0*rxT);         % weight x receive part of the remodulation
txTime = @(x, z, th) (x*sin(th) + z*cos(th))/c;

% f-k jobs: [map (1 ERM, 2 exact), FFT length, nearest (1) or linear (0), time-centering]
jobMain = [1 ntFFT 0 1; 2 ntFFT 0 1];
jobTest = [2 2048 1 0; 2 2048 0 0; 2 4096 0 0; 2 4096 0 1; 2 8192 0 1];

%% Simulate, then beamform every angle three ways
img  = zeros(numel(zg), numel(xg), Nang, 3, 2);     % (z, x, angle, method, phantom)
imgT = zeros(numel(zg), numel(xg), size(jobTest,1), 2);
tDAS = 0;  tFK = zeros(1,2);
for ia = 1:Nang
    th = theta(ia);
    S = cell(1,2);
    for ph = 1:2                                 % channel data (IQ, nt x Nel)
        if ph == 1, X0 = ptsX; Z0 = ptsZ; A0 = ptsA; else, X0 = spX; Z0 = spZ; A0 = spA; end
        tau = txTime(X0, Z0, th) + sqrt((X0 - xe).^2 + Z0.^2)/c;   % two-way delay (us)
        w   = A0 .* exp(-1i*2*pi*f0*tau);        % exact carrier phase
        q   = tau*fs;  i0 = floor(q);  fr = q - i0;
        li  = i0 + 1 + nt*(0:Nel-1);
        Sd  = accumarray(li(:), w(:).*(1-fr(:)), [nt*Nel 1]) + accumarray(li(:)+1, w(:).*fr(:), [nt*Nel 1]);
        S{ph} = conv2(reshape(Sd, nt, Nel), envK, 'same');
    end
    % ---- DAS
    t1 = tic;
    T  = txTime(Px(pp), Pz(pp), th) + rxT;  q = T*fs;  i0 = floor(q);  fr = q - i0;
    li = i0 + 1 + nt*(ee - 1);  txRot = exp(1i*2*pi*f0*txTime(Px, Pz, th));
    for ph = 1:2
        v = ((1-fr).*S{ph}(li) + fr.*S{ph}(li+1)) .* Wc;
        img(:,:,ia,1,ph) = reshape(accumarray(pp, v, [nPix 1]) .* txRot, numel(zg), numel(xg));
    end
    tDAS = tDAS + toc(t1);
    % ---- f-k
    jobs = jobMain;  if ia == iMid, jobs = [jobMain; jobTest]; end
    for jb = 1:size(jobs, 1)
        mode = jobs(jb,1);  ntq = jobs(jb,2);  nearest = jobs(jb,3);  tc = jobs(jb,4)*tCenter;
        t1 = tic;
        Dfq  = fs/ntq;  fq = [0:ntq/2-1, -ntq/2:-1]*Dfq + f0;     % RF frequency of each bin
        bq   = find(fq >= 0.8 & fq <= 9.2);  [fB, o] = sort(fq(bq));  bq = bq(o);
        fB   = fB(:);  nB = numel(fB);
        if mode == 1        % Garcia ERM fit
            v  = c/sqrt(1 + cos(th) + sin(th)^2);  be = (1 + cos(th))^1.5/(1 + cos(th) + sin(th)^2);
            ga = sin(th)/(2 - cos(th));  kzh = (KZ - ga*KX)/be;
            F  = v*sqrt(KX.^2 + kzh.^2);  J = v*kzh./(be*sqrt(KX.^2 + kzh.^2) + eps);  ok = kzh > 0;
        else                % exact transmit + receive wavevector map
            M  = KX*sin(th) + KZ*cos(th);  N2 = KX.^2 + KZ.^2;
            F  = c*N2./(2*M + eps);  J = abs(c*(2*KZ.*M - N2*cos(th))./(2*M.^2 + eps));
            ok = M > 0 & (KZ - F*cos(th)/c) >= 0;
        end
        q  = (F - fB(1))/Dfq + 1;
        if nearest, i0 = round(q); w = zeros(size(q)); else, i0 = floor(q); w = q - i0; end
        ok = ok & i0 >= 1 & i0 < nB;
        li = i0(ok) + nB*(ceil(find(ok)/NZ) - 1);  w = w(ok);
        Wt = J(ok) .* exp(-1i*2*pi*F(ok)*tc);    % Jacobian, undo the time-centering
        for ph = 1:2
            D  = fft(S{ph}, ntq, 1);
            D  = D(bq, :) .* exp(1i*2*pi*fB*tc);         % RF frequency rows, time-centered
            D  = fft(D, nxp, 2);                         % data wavenumber kx_d
            if useFan, D = D .* (abs(kxd) <= fB*sin(alphaMax)/c); end
            Du = zeros(nB, nxf);  Du(:, 1:nxp/2) = D(:, 1:nxp/2);  Du(:, end-nxp/2+1:end) = D(:, nxp/2+1:end);
            d  = ifft(Du, [], 2) .* exp(1i*2*pi*(fB*sin(th)/c)*xf);   % steering compensation
            D2 = fft(d, [], 2);                          % nB x nxf, image Kx
            I  = zeros(NZ, nxf);
            I(ok) = ((1 - w).*D2(li) + w.*D2(li+1)) .* Wt;
            im = ifft2(I);
            if jb <= 2, img(:,:,ia,mode+1,ph) = im(iz, ix); else, imgT(:,:,jb-2,ph) = im(iz, ix); end
        end
        if jb <= 2, tFK(mode) = tFK(mode) + toc(t1); end
    end
end
fprintf('time per angle (2 phantoms): DAS %.2f s, f-k ERM %.2f s, f-k exact %.2f s -> f-k %.1fx faster\n', ...
    tDAS/Nang, tFK(1)/Nang, tFK(2)/Nang, tDAS/mean(tFK));
fprintf('work per angle and phantom: DAS %d (pixel, element) pairs; f-k FFTs %dx%d + 2 x %dx%d, %d interpolations\n\n', ...
    numel(pp), ntFFT, Nel, nB, nxf, nnz(ok));

%% Normalize each method by its peak on the central point target (0 deg)
names = {'DAS', 'f-k ERM', 'f-k exact'};
nrm = @(I3) max(max(abs(I3(abs(zg-15) < 1, abs(xg) < 1))));
for m = 1:3, sc = nrm(img(:,:,iMid,m,1));  img(:,:,:,m,:) = img(:,:,:,m,:)/sc; end
for k = 1:size(jobTest,1), sc = nrm(imgT(:,:,k,1));  imgT(:,:,k,:) = imgT(:,:,k,:)/sc; end

%% Metrics on the central target (0, 15 mm)
xl = -1.5:0.005:1.5;  zl = -1.2:0.002:1.2;
box = abs(XG) <= 4 & abs(ZG - 15) <= 4;  mlobe = abs(XG) <= 1.0 & abs(ZG - 15) <= 0.45;
sets = {iMid, 1:Nang, find(thetaDeg == 9)};  setName = {'0 deg', sprintf('%d angles', Nang), '9 deg'};
fprintf('%-10s %-9s %8s %8s %7s %9s %9s %8s\n', 'method', 'set', 'lat-6dB', 'ax-6dB', 'PSL', 'dx peak', 'dz peak', 'A(22)/');
fprintf('%-10s %-9s %8s %8s %7s %9s %9s %8s\n', '', '', '(mm)', '(mm)', '(dB)', '(um)', '(um)', 'A(15)');
for s = 1:3
    for m = 1:3
        E  = abs(sum(img(:,:,sets{s},m,1), 3));
        sub = E .* box;  [pk, k] = max(sub(:));  [kz, kx] = ind2sub(size(E), k);
        rz = max(1,kz-12):min(numel(zg),kz+12);  rx = max(1,kx-8):min(numel(xg),kx+8);
        Lc = interp1(zg(rz), E(rz, kx), zg(kz) + zl, 'spline');  [~, b] = max(Lc);  zpk = zg(kz) + zl(b);
        rowv = interp1(zg(rz), E(rz, rx), zpk, 'spline');
        Lr = interp1(xg(rx), rowv, xg(kx) + xl, 'spline');  [pr, a] = max(Lr);  xpk = xg(kx) + xl(a);
        colv = interp1(xg(rx), E(rz, rx).', xpk, 'spline');
        Lc = interp1(zg(rz), colv(:), zg(kz) + zl, 'spline');  [pc, b] = max(Lc);  zpk = zg(kz) + zl(b);
        wl = sum(Lr >= pr/2)*0.005;  wa = sum(Lc >= pc/2)*0.002;
        PSL = 20*log10(max(E(box & ~mlobe))/pk);
        A22 = max(max(E(abs(ZG-22) < 1 & abs(XG) < 1)))/pk;
        fprintf('%-10s %-9s %8.3f %8.3f %7.1f %9.1f %9.1f %8.2f\n', names{m}, setName{s}, wl, wa, PSL, xpk*1e3, (zpk-15)*1e3, A22);
    end
end
fprintf('theory: receive-only 1.207 lam F# = %.3f mm; axial -6 dB of the envelope 2.355 sig_t c/2 = %.3f mm\n', ...
    1.207*lam*Fnum, 2*sqrt(2*log(2))*sig_t*c/2);
fprintf('pixel grid: dx %.0f um, dz %.0f um (peak positions are spline-refined)\n\n', dxf*1e3, (zg(2)-zg(1))*1e3);

%% Speckle: cyst contrast, clutter beside the block, residual against DAS
rr   = sqrt((XG-cx).^2 + (ZG-cz).^2);
cin  = rr <= R - 0.5;  cout = rr >= R + 1 & rr <= R + 3.5 & ZG <= 38.5;
blk  = abs(XG) < 7 & ZG > 28 & ZG < 38 & rr > R + 1;   side = XG > 9.5 & ZG > 29 & ZG < 38;
gres = @(I, ref) 20*log10(norm((sum(conj(I).*ref, 2)./sum(abs(I).^2, 2)).*I - ref, 'fro')/norm(ref, 'fro'));
spk  = @(I) [10*log10(mean(abs(I(cin)).^2)/mean(abs(I(cout)).^2)), 10*log10(mean(abs(I(side)).^2)/mean(abs(I(blk)).^2))];
fprintf('%-10s %-9s %10s %10s %12s\n', 'method', 'set', 'contrast', 'beside', 'residual vs');
fprintf('%-10s %-9s %10s %10s %12s\n', '', '', '(dB)', 'block (dB)', 'DAS (dB)');
for s = 1:2
    ref = sum(img(:,:,sets{s},1,1) + img(:,:,sets{s},1,2), 3);
    for m = 1:3
        Ic = sum(img(:,:,sets{s},m,1) + img(:,:,sets{s},m,2), 3);  r = spk(sum(img(:,:,sets{s},m,2), 3));
        if m == 1, rs = '(reference)'; else, rs = sprintf('%.1f', gres(Ic, ref)); end
        fprintf('%-10s %-9s %10.1f %10.1f %12s\n', names{m}, setName{s}, r(1), r(2), rs);
    end
end
fprintf('(residual after one complex gain per depth row: removes the depth-gain difference\n');
fprintf(' between unit-sum DAS and Stolt; what remains is the k-space weighting difference)\n\n');

%% Interpolation along f, zero-padding and time-centering (0 deg, exact map)
fprintf('%-8s %6s %8s %9s %10s %10s %12s\n', 'interp', 'ntFFT', 'df(kHz)', 'centered', 'contrast', 'beside', 'residual vs');
fprintf('%-8s %6s %8s %9s %10s %10s %12s\n', '', '', '', '', '(dB)', 'block (dB)', 'DAS (dB)');
ref = img(:,:,iMid,1,1) + img(:,:,iMid,1,2);
for k = 1:size(jobTest,1)
    Ic = imgT(:,:,k,1) + imgT(:,:,k,2);  r = spk(imgT(:,:,k,2));
    nm = 'linear'; if jobTest(k,3), nm = 'nearest'; end
    yn = 'no'; if jobTest(k,4), yn = 'yes'; end
    fprintf('%-8s %6d %8.2f %9s %10.1f %10.1f %12.1f\n', nm, jobTest(k,2), fs/jobTest(k,2)*1e3, yn, r(1), r(2), gres(Ic, ref));
end
fprintf('\ntotal run time %.1f s\n', toc(tStart));

%% Figures
dr = 50;  dB = @(I) max(20*log10(abs(I)/max(abs(I(:))) + eps), -dr);
figure('Position', [40 40 1150 720]);
for m = 1:3
    Ic = sum(img(:,:,:,m,1) + img(:,:,:,m,2), 3);
    subplot(2,3,m); imagesc(xg, zg, dB(Ic), [-dr 0]); axis image; colormap(gray);
    title(sprintf('%s, %d angles (dB)', names{m}, Nang)); xlabel('x (mm)'); ylabel('z (mm)');
end
subplot(2,3,4);
E1 = abs(img(:,:,iMid,1,1));  E2 = abs(img(:,:,iMid,3,1));  E3 = abs(sum(img(:,:,:,3,1),3));
[~, jz] = min(abs(zg - 15));
plot(xg, 20*log10(E1(jz,:)/max(E1(jz,:))), 'k', xg, 20*log10(E2(jz,:)/max(E2(jz,:))), 'r--', ...
     xg, 20*log10(E3(jz,:)/max(E3(jz,:))), 'b');
axis([-12 12 -60 0]); grid on; xlabel('x (mm)'); ylabel('dB');
legend('DAS, 0 deg', 'f-k exact, 0 deg', 'f-k exact, 7 angles'); title('lateral profile, z = 15 mm');
subplot(2,3,5);
Ic = imgT(:,:,2,1) + imgT(:,:,2,2);
imagesc(xg, zg, dB(Ic), [-dr 0]); axis image; title('0 deg, linear, 2048 pts, no centering');
xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,6);
Ic = imgT(:,:,4,1) + imgT(:,:,4,2);
imagesc(xg, zg, dB(Ic), [-dr 0]); axis image; title('0 deg, linear, 4096 pts, centered');
xlabel('x (mm)'); ylabel('z (mm)');

What to look for: widths 0.535 / 0.565 / 0.555 mm at 0° (DAS / ERM / exact), all 0.505 mm with 7 angles; PSL −35.1 / −33.6 / −34.2 dB compounded; at 9° the ERM peak sits 35 µm off, the exact map on target; cyst contrast −26.5 / −26.9 / −27.8 dB; the interpolation table drops from −10.6 dB (2048, uncentred) to −18.2 dB beside the block (DAS −18.7 dB).

m24_fk_migration_zero_offset.mExploding reflector: Stolt with v = c/2 against SAFT, v = c failure, the λ/4 pitch rule · change pitches, Fnum, targets · runtime ≈ 10 s

%% m24_fk_migration_zero_offset.m  --  exploding reflector and Stolt migration of zero-offset data
% Purpose : zero-offset (monostatic synthetic aperture) data: every element fires alone and
%           records only its own echo, so the two-way delay to a point is 2r/c. That is the
%           one-way delay of a source at the point radiating at speed c/2 (the exploding
%           reflector). Stolt migration with v = c/2 maps the data spectrum along
%           f = (c/2) sqrt(Kx^2 + Kz^2) with Jacobian (c/2) Kz/sqrt(Kx^2 + Kz^2).
%           The script checks (1) Stolt against delay-and-sum (SAFT) on the same pixels,
%           (2) what using v = c instead of c/2 does, (3) the lambda/4 pitch rule for
%           zero-offset data: the data occupy |kx| <= 2f/c, so the pitch must satisfy
%           p <= c/(4 f_max) to avoid spatial aliasing at all angles.
% Change  : pitch list, v factor, target positions, f0.
% Expect  : Stolt and SAFT peaks at the true positions (errors of a few um), the same
%           lateral width (two-way, full 38.4 mm aperture: ~0.1-0.2 mm) within a few %.
%           v = c puts the point at twice its depth and defocuses it. Artifact level away
%           from the targets drops sharply from pitch 0.30 mm (lambda) to 0.075 mm (lambda/4).
% Model   : linear acoustics, Born, point-like elements, no attenuation or spreading,
%           Gaussian pulse, IQ at 20 MHz. Runtime ~10 s in Octave.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);
tStart = tic;

%% Parameters
c = 1.540; f0 = 5; lam = c/f0; FBW = 0.6; fs = 20;   % mm/us, MHz, mm, -, MHz
L = 38.4;                                            % aperture (mm), same for every pitch
pitches = [0.30 0.15 0.075];                         % lambda, lambda/2, lambda/4 at 5 MHz
ptsX = [0 -6 6 0].';  ptsZ = [10 20 20 30].';
nt = 1100;  ntFFT = 4096;  tc = nt/(2*fs);           % samples, FFT length, centering (us)
sig_f = FBW*f0/(2*sqrt(2*log(2)));  sig_t = 1/(2*pi*sig_f);
kt = (-ceil(4*sig_t*fs):ceil(4*sig_t*fs)).';  envK = exp(-(kt/fs).^2/(2*sig_t^2));
dxf = 0.075;                                         % common lateral image grid (mm)
Fnum = 1.5;  alphaMax = atan(1/(2*Fnum));           % receive F# (SAFT aperture, Stolt fan)
NZ = 1024;  dKz = 2*(fs/2048)/c;  Kz = (0:NZ-1).'*dKz;  zf = (0:NZ-1).'/(NZ*dKz);
Dfq = fs/ntFFT;  fq = [0:ntFFT/2-1, -ntFFT/2:-1]*Dfq + f0;
bq = find(fq >= 0.8 & fq <= 9.2);  [fB, o] = sort(fq(bq));  bq = bq(o);  fB = fB(:);  nB = numel(fB);
fprintf('zero-offset aliasing-free pitch c/(4 f) = %.3f mm at f0, %.3f mm at %.1f MHz (-24 dB band edge)\n\n', ...
    c/(4*f0), c/(4*8), 8);

fprintf('%-6s %-6s %-6s %8s %8s %9s %9s %10s\n', 'pitch', 'p/lam', 'method', 'target', 'lat-6dB', 'dx peak', 'dz peak', 'artifact');
fprintf('%-6s %-6s %-6s %8s %8s %9s %9s %10s\n', '(mm)', '', '', '(x,z)', '(mm)', '(um)', '(um)', '(dB)');
store = cell(1, numel(pitches));
for ip = 1:numel(pitches)
    p = pitches(ip);  Nel = round(L/p);  xe = ((0:Nel-1) - (Nel-1)/2)*p;
    %% zero-offset channel data, IQ (nt x Nel)
    tau = 2*sqrt((ptsX - xe).^2 + ptsZ.^2)/c;            % two-way delay (us)
    w = exp(-1i*2*pi*f0*tau);  q = tau*fs;  i0 = floor(q);  fr = q - i0;
    li = i0 + 1 + nt*(0:Nel-1);
    Sd = accumarray(li(:), w(:).*(1-fr(:)), [nt*Nel 1]) + accumarray(li(:)+1, w(:).*fr(:), [nt*Nel 1]);
    S  = conv2(reshape(Sd, nt, Nel), envK, 'same');
    %% Stolt (v = c/2), lateral zero-padding to 4x the aperture, upsampled to dxf
    nxp = 2^nextpow2(4*Nel);  U = round(p/dxf);  nxf = U*nxp;
    xf  = xe(1) + (0:nxf-1)*dxf;  Kx = [0:nxf/2-1, -nxf/2:-1]/(nxf*dxf);
    [KX, KZ] = meshgrid(Kx, Kz);
    D  = fft(S, ntFFT, 1);  D = D(bq, :) .* exp(1i*2*pi*fB*tc);  D = fft(D, nxp, 2);
    kxd = [0:nxp/2-1, -nxp/2:-1]/(nxp*p);  Dpre = D;
    D  = D .* (abs(kxd) <= 2*fB*sin(alphaMax)/c);       % fan filter: zero-offset kx = 2 f sin(alpha)/c
    Du = zeros(nB, nxf);  Du(:, 1:nxp/2) = D(:, 1:nxp/2);  Du(:, end-nxp/2+1:end) = D(:, nxp/2+1:end);
    imgs = cell(1,2);
    for vf = [0.5 1]                                     % velocity c/2 (correct) and c (wrong)
        v = vf*c;  F = v*sqrt(KX.^2 + KZ.^2);  J = v*KZ./(sqrt(KX.^2 + KZ.^2) + eps);
        qq = (F - fB(1))/Dfq + 1;  k0 = floor(qq);  ww = qq - k0;  ok = k0 >= 1 & k0 < nB;
        lk = k0(ok) + nB*(ceil(find(ok)/NZ) - 1);
        I = zeros(NZ, nxf);  I(ok) = ((1 - ww(ok)).*Du(lk) + ww(ok).*Du(lk+1)) .* J(ok) .* exp(-1i*2*pi*F(ok)*tc);
        imgs{1 + (vf == 1)} = ifft2(I);
    end
    E = abs(imgs{1});  E = E/max(E(:));
    store{ip} = struct('E', E, 'xf', xf, 'Ew', abs(imgs{2})/max(abs(imgs{1}(:))), 'S', S, 'xe', xe, 'D', Dpre);
    %% SAFT (delay and sum, dynamic aperture F#) on the Stolt pixels around each target
    for it = 1:numel(ptsX)
        jx = find(abs(xf - ptsX(it)) <= 3);  jz = find(abs(zf - ptsZ(it)) <= 1);
        [XL, ZL] = meshgrid(xf(jx), zf(jz));
        Td = 2*sqrt((XL(:) - xe).^2 + ZL(:).^2)/c;       % pixels x elements (us)
        qd = Td*fs;  kd = floor(qd);  fd = qd - kd;
        okd = kd >= 0 & kd < nt - 1 & abs(XL(:) - xe) <= ZL(:)/(2*Fnum);
        kd(~okd) = 0;  ld = kd + 1 + nt*(0:Nel-1);
        V  = ((1-fd).*S(ld) + fd.*S(ld+1)) .* okd .* exp(1i*2*pi*f0*Td);
        Bd = reshape(abs(sum(V, 2)), size(XL));
        Bs = E(jz, jx);
        for m = 1:2
            if m == 1, B = Bs; nm = 'Stolt'; else, B = Bd; nm = 'SAFT'; end
            B = B/max(B(:));  [~, k] = max(B(:));  [kz, kx] = ind2sub(size(B), k);
            xs = XL(1,:);  zs = ZL(:,1);  xq = xs(kx) + (-0.6:0.002:0.6);  zq = zs(kz) + (-0.4:0.002:0.4);
            Lc = interp1(zs, B(:, kx), zq, 'spline');  [~, b2] = max(Lc);  zpk = zq(b2);
            rowv = interp1(zs, B, zpk, 'spline');  Lr = interp1(xs, rowv, xq, 'spline');
            [pr, a2] = max(Lr);  xpk = xq(a2);  wl = sum(Lr >= pr/2)*0.002;
            far = abs(XL - xpk) > 1.0 | abs(ZL - zpk) > 0.5;
            fprintf('%-6.3f %-6.2f %-6s %8s %8.3f %9.1f %9.1f %10.1f\n', p, p/lam, nm, sprintf('%g,%g', ptsX(it), ptsZ(it)), ...
                wl, (xpk - ptsX(it))*1e3, (zpk - ptsZ(it))*1e3, 20*log10(max(B(far))));
        end
    end
    %% wrong velocity: where does the target at 10 mm go?
    Ew = store{ip}.Ew;  sel = abs(xf) < 3;  [pw, k] = max(max(Ew(zf < 26, sel), [], 2));
    fprintf('       v = c instead of c/2: brightest point near x = 0 at z = %.2f mm (true 10), peak %.1f dB re correct\n\n', zf(k), 20*log10(pw));
end
fprintf('uniform kx support |kx| <= 2 sin(alpha)/lam gives a -6 dB width 1.207 lam/(4 sin alpha) = %.3f mm\n', 1.207*lam/(4*sin(alphaMax)));
fprintf('fan edge reaches the pitch Nyquist 1/(2p) at f = c/(4 p sin alpha): %.1f MHz (p = 0.30), %.1f MHz (0.15), %.1f MHz (0.075)\n', ...
    c./(4*pitches*sin(alphaMax)));
fprintf('total run time %.1f s\n', toc(tStart));

%% Figures
s = store{1};  t = (0:nt-1)/fs;
figure('Position', [50 50 1150 650]);
subplot(2,3,1); imagesc(s.xe, t, real(s.S)); colormap(gray); xlabel('element x (mm)'); ylabel('t (us)');
title('zero-offset data (real IQ), pitch 0.30 mm');
subplot(2,3,2); nxp = size(s.D,2); kxd = [0:nxp/2-1, -nxp/2:-1]/(nxp*pitches(1));
P = 20*log10(abs(fftshift(s.D, 2)) + eps);  P = P - max(P(:));
imagesc(fftshift(kxd), fB, max(P, -50)); axis xy; hold on; axis([-1.8 1.8 fB(1) fB(end)]);
plot(2*fB/c, fB, 'r--', -2*fB/c, fB, 'r--'); plot([1 1]/(2*pitches(1)), fB([1 end]), 'y', -[1 1]/(2*pitches(1)), fB([1 end]), 'y');
xlabel('k_x (cycles/mm)'); ylabel('f (MHz)'); title('|D(k_x,f)| (dB), p = 0.30 mm: red |k_x| = 2f/c (folds), yellow Nyquist');
subplot(2,3,3); s = store{1}; imagesc(s.xf, zf, max(20*log10(s.E + eps), -50), [-50 0]); axis image; axis([-10 10 5 35]);
title('Stolt, v = c/2, pitch 0.30 mm (dB)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,4); s = store{3}; imagesc(s.xf, zf, max(20*log10(s.E + eps), -50), [-50 0]); axis image; axis([-10 10 5 35]);
title('Stolt, v = c/2, pitch 0.075 mm (dB)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,5); s = store{3}; imagesc(s.xf, zf, max(20*log10(s.Ew/max(s.Ew(:)) + eps), -50), [-50 0]); axis image; axis([-15 15 5 70]);
title('Stolt with v = c (wrong): depth x2, defocused'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,6); s1 = store{1}; s3 = store{3}; [~, jz] = min(abs(zf - 20));
plot(s1.xf, 20*log10(s1.E(jz,:) + eps), 'r', s3.xf, 20*log10(s3.E(jz,:) + eps), 'k'); axis([-12 12 -60 0]); grid on;
xlabel('x (mm)'); ylabel('dB'); legend('pitch 0.30 mm', 'pitch 0.075 mm'); title('lateral profile at z = 20 mm');

What to look for: at 0.15 and 0.075 mm pitch Stolt and SAFT agree (0.28 against 0.29 mm width, uniform-\(k_x\) theory 0.294 mm; positions within 2 µm); at 0.30 mm the fan is clipped by Nyquist above 4.1 MHz, Stolt widens to 0.38 mm and aliasing artifacts rise to −21.5 dB; \(v = c\) images the 10 mm target at 19.94 mm.

m24_kspace_evanescent.mThe data cone, energy beyond Nyquist versus depth, evanescent gains, phase-shift depth loop against Stolt · change depths, pitch · runtime ≈ 3 s

%% m24_kspace_evanescent.m  --  channel data in (kx, f): the propagating cone, aliasing, evanescent waves
% Purpose : (1) 2D FFT of 0-degree plane-wave channel data. The receive leg is a one-way
%           wave field on the array, so its energy lies inside the cone |kx| <= f/c; the
%           region outside is evanescent and empty. (2) Spatial aliasing: with pitch p the
%           kx axis wraps at 1/(2p); the energy that would fold is measured per target
%           depth for p = 0.30 mm (about lambda at 5 MHz). (3) Why evanescent components are
%           never back-propagated: their gain exp(2 pi z sqrt(kx^2 - (f/c)^2)) explodes.
%           (4) Phase-shift imaging (depth loop, exact plane-wave map, no interpolation)
%           against the Stolt version of the same map (one interpolation): equal images.
% Change  : target depths, pitch, the depth step of the phase-shift loop.
% Expect  : energy outside the cone below -40 dB (window leakage only); the share of
%           energy beyond |kx| = 1/(2 x 0.30 mm) grows from a few % for deep targets to
%           ~30 % for a target at 5 mm; Stolt vs phase-shift residual near -30 dB.
% Model   : as m24_stolt_planewave.m (Born, point elements, no attenuation), pitch 0.15 mm
%           (lambda/2 at 5 MHz, 256 elements) so the reference spectrum is not aliased.
%           Runtime ~10 s in Octave.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);
tStart = tic;

%% Parameters and data (0-degree plane wave, IQ)
c = 1.540; f0 = 5; lam = c/f0; FBW = 0.6; fs = 20;
p = 0.15;  Nel = 256;  xe = ((0:Nel-1) - (Nel-1)/2)*p;        % 38.4 mm aperture
depths = [5 10 20 30];  nt = 1000;
sig_f = FBW*f0/(2*sqrt(2*log(2)));  sig_t = 1/(2*pi*sig_f);
kt = (-ceil(4*sig_t*fs):ceil(4*sig_t*fs)).';  envK = exp(-(kt/fs).^2/(2*sig_t^2));
nfft = 1024;  Df = fs/nfft;  fq = [0:nfft/2-1, -nfft/2:-1]*Df + f0;
bq = find(fq >= 0.8 & fq <= 9.2);  [fB, o] = sort(fq(bq));  bq = bq(o);  fB = fB(:);  nB = numel(fB);
nxp = 512;  kx = [0:nxp/2-1, -nxp/2:-1]/(nxp*p);              % cycles/mm
spec = cell(1, numel(depths));  Sall = 0;
tw = 0.5 - 0.5*cos(pi*min(1, min(0:Nel-1, Nel-1:-1:0)/(0.125*(Nel-1))));   % Tukey(0.25) taper
fprintf('%8s %14s %14s %16s %18s\n', 'depth', 'outside cone', 'outside cone', 'beyond 1/(2x0.3)', 'max receive angle');
fprintf('%8s %14s %14s %16s %18s\n', '(mm)', 'rect (dB)', 'Tukey (dB)', '(% of energy)', 'to array edge');
for id = 1:numel(depths)
    z0 = depths(id);  tau = (z0 + sqrt(xe.^2 + z0^2))/c;      % 0-degree plane wave, two-way
    q = tau*fs;  i0 = floor(q);  fr = q - i0;  w = exp(-1i*2*pi*f0*tau);
    Sd = zeros(nt, Nel);
    Sd(sub2ind([nt Nel], i0 + 1, 1:Nel)) = w.*(1-fr);  Sd(sub2ind([nt Nel], i0 + 2, 1:Nel)) = Sd(sub2ind([nt Nel], i0 + 2, 1:Nel)) + w.*fr;
    S  = conv2(Sd, envK, 'same');  Sall = Sall + S;
    D  = fft(S, nfft, 1);  D = fft(D(bq, :), nxp, 2);  P = abs(D).^2;
    Dw = fft(S, nfft, 1);  Dw = fft(Dw(bq, :).*tw, nxp, 2);  Pw = abs(Dw).^2;   % tapered aperture
    out  = abs(kx) > fB/c;  wrap = abs(kx) > 1/(2*0.30);
    spec{id} = D;
    fprintf('%8.0f %14.1f %14.1f %16.1f %16.1f deg\n', z0, 10*log10(sum(P(out))/sum(P(:))), 10*log10(sum(Pw(out))/sum(Pw(:))), ...
        100*sum(P(wrap & ~out))/sum(P(:)), atan(max(abs(xe))/z0)*180/pi);
end
fprintf('kx Nyquist: %.2f /mm (p = 0.15 mm), %.2f /mm (p = 0.30 mm); cone edge f/c at 5 MHz: %.2f /mm\n', ...
    1/(2*p), 1/(2*0.30), f0/c);
fprintf('the cone reaches the p = 0.30 mm Nyquist at f = c/(2p) = %.2f MHz: above it, steep arrivals alias\n\n', c/(2*0.30));

%% Evanescent gain if one tried to back-propagate
fprintf('back-propagation gain of an evanescent component at 5 MHz, exp(2 pi z sqrt(kx^2-(f/c)^2)):\n');
fprintf('%12s %12s %12s %12s\n', 'kx/(f/c)', 'z = lam/4', 'z = 1 mm', 'z = 10 mm');
for r = [1.01 1.1 1.5]
    a = 2*pi*sqrt(r^2 - 1)*f0/c;                               % decay rate (1/mm)
    fprintf('%12.2f %10.1f dB %10.1f dB %10.0f dB\n', r, 20*log10(exp(a*lam/4)), 20*log10(exp(a*1)), 20*log10(exp(a*10)));
end
fprintf('(the same numbers, negative, are how fast an evanescent field decays away from its source)\n\n');

%% Phase-shift imaging vs Stolt, both with the exact 0-degree map Kz = kz + f/c
D = fft(Sall, nfft, 1);  D = fft(D(bq, :), nxp, 2);            % nB x nxp
dz = 0.077;  zg = (round(3/dz):round(33/dz)).'*dz;  Nz = numel(zg);   % on the Stolt depth grid
KZr = sqrt(max((fB/c).^2 - kx.^2, 0));  prop = (fB/c).^2 >= kx.^2;  % receive kz, propagating mask
t1 = tic;
psi = D .* prop .* exp(1i*2*pi*(KZr + fB/c)*zg(1));            % field at the first depth
step = exp(1i*2*pi*(KZr + fB/c)*dz) .* prop;
imPS = zeros(Nz, nxp);
for k = 1:Nz
    imPS(k, :) = sum(psi, 1);                                  % imaging condition t = 0: sum over f
    psi = psi .* step;                                         % one depth step, transmit + receive
end
imPS = ifft(imPS, [], 2);                                      % kx -> x
tPS = toc(t1);
% Stolt with the same map (interpolation along f instead of the depth loop)
t1 = tic;
ntq = 4096;  Dq = fs/ntq;  fq2 = [0:ntq/2-1, -ntq/2:-1]*Dq + f0;
b2 = find(fq2 >= 0.8 & fq2 <= 9.2);  [fB2, o] = sort(fq2(b2));  b2 = b2(o);  fB2 = fB2(:);  n2 = numel(fB2);
tc = nt/(2*fs);
D2 = fft(Sall, ntq, 1);  D2 = fft(D2(b2, :) .* exp(1i*2*pi*fB2*tc), nxp, 2);
NZ = 2048;  dKz = 1/(NZ*dz);  Kz = (0:NZ-1).'*dKz;  zS = (0:NZ-1).'*dz;
[KX, KZ] = meshgrid(kx, Kz);
F = c*(KX.^2 + KZ.^2)./(2*KZ + eps);  J = abs(c*(KZ.^2 - KX.^2)./(2*KZ.^2 + eps));
ok = KZ > 0 & (KZ - F/c) >= 0;
qq = (F - fB2(1))/Dq + 1;  k0 = floor(qq);  ww = qq - k0;  ok = ok & k0 >= 1 & k0 < n2;
lk = k0(ok) + n2*(ceil(find(ok)/NZ) - 1);
I = zeros(NZ, nxp);  I(ok) = ((1 - ww(ok)).*D2(lk) + ww(ok).*D2(lk+1)) .* J(ok) .* exp(-1i*2*pi*F(ok)*tc);
imS = ifft2(I);  [~, j0] = min(abs(zS - zg(1)));  imS = imS(j0:j0+Nz-1, :);
tS = toc(t1);
% The phase-shift sum over f needs no weight: Stolt's Jacobian df = |df/dKz| dKz is exactly
% what turns the same integral over f into an integral over Kz.
fprintf('phase-shift depth loop: %.2f s (%d depths); Stolt: %.2f s\n', tPS, Nz, tS);
sel = abs(((0:nxp-1)*p + xe(1))) <= 15;  A = imPS(:, sel);  B = imS(:, sel);
al = (A(:)'*B(:))/(A(:)'*A(:));
E1 = abs(al*A);  E2 = abs(B);  m = E2 > 0.01*max(E2(:));
xs = (0:nxp-1)*p + xe(1);  [~, kx0] = min(abs(xs));
fprintf('Stolt (with Jacobian) vs phase-shift (no weight), same pixels:\n');
for z0 = depths
    [~, kz] = min(abs(zg - z0));  a1 = al*imPS(kz, kx0);  b1 = imS(kz, kx0);
    fprintf('  target at %2g mm: amplitude ratio %.3f, phase difference %.3f rad\n', z0, abs(b1)/abs(a1), angle(b1/a1));
end
fprintf('  max |difference| %.1f dB re peak; mean |dB difference| over pixels above -40 dB: %.2f dB\n', ...
    20*log10(max(abs(al*A(:) - B(:)))/max(E2(:))), mean(abs(20*log10(E1(m)./E2(m)))));
fprintf('total run time %.1f s\n', toc(tStart));

%% Figures
figure('Position', [50 50 1100 650]);
for id = [1 3]
    subplot(2,3,1 + (id == 3));
    P = 20*log10(abs(fftshift(spec{id}, 2)) + eps);  P = P - max(P(:));
    imagesc(fftshift(kx), fB, max(P, -50), [-50 0]); axis xy; hold on;
    plot(fB/c, fB, 'r--', -fB/c, fB, 'r--');
    plot([1 1]/(2*0.30), fB([1 end]), 'y:', -[1 1]/(2*0.30), fB([1 end]), 'y:');
    xlabel('k_x (cycles/mm)'); ylabel('f (MHz)'); title(sprintf('|D(k_x,f)| dB, target at %g mm', depths(id)));
end
colormap(gray);
subplot(2,3,3);
kr = linspace(1, 1.5, 200);  plot(kr, -20*log10(exp(2*pi*sqrt(kr.^2 - 1)*f0/c*lam)), 'k'); grid on;
xlabel('k_x / (f/c)'); ylabel('dB per wavelength'); title('evanescent decay over one wavelength, 5 MHz');
xs = (0:nxp-1)*p + xe(1);
subplot(2,3,4); imagesc(xs, zg, max(20*log10(abs(imPS)/max(abs(imPS(:)))), -50), [-50 0]); axis image; axis([-15 15 3 33]);
title('phase-shift (depth loop), dB'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,5); imagesc(xs, zg, max(20*log10(abs(imS)/max(abs(imS(:)))), -50), [-50 0]); axis image; axis([-15 15 3 33]);
title('Stolt (one interpolation), dB'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,6); imagesc(xs, zg, max(20*log10(abs(al*imPS - imS)/max(abs(imS(:)))), -60), [-60 0]); axis image; axis([-15 15 3 33]);
colorbar; title('|difference| (dB re peak)'); xlabel('x (mm)'); ylabel('z (mm)');

What to look for: energy outside the cone −22 to −35 dB with a rectangular aperture, −29 to −53 dB with a Tukey taper (leakage, not evanescent waves); 82 / 65 / 34 / 10 % of the energy beyond the 0.30 mm Nyquist for targets at 5 / 10 / 20 / 30 mm; Stolt and phase-shift peaks within 2.1 % in amplitude and 0.000 rad in phase.

m24_phase_shift_layered.m8 mm fat-like layer: phase-shift migration against single-speed Stolt, ray-traced data · change c1, d1, targets · runtime ≈ 3 s

%% m24_phase_shift_layered.m  --  phase-shift migration through a layer vs Stolt with one speed
% Purpose : a 0-degree plane wave images point targets below an 8 mm fat-like layer
%           (c1 = 1450 m/s) on tissue (c2 = 1540 m/s). Channel data use ray-traced (Fermat)
%           receive delays and the exact vertical transmit delay. Three reconstructions:
%             phase-shift : depth loop with the local speed, step exp(j2 pi (kz(z) + f/c(z)) dz),
%                           kz(z) = sqrt((f/c(z))^2 - kx^2); kx is conserved across a flat
%                           interface (Snell), so this is exact for horizontal layers
%             Stolt 1540  : exact 0-degree map with c = 1540 m/s everywhere
%             Stolt cbar  : the same with the time-average speed down to 25 mm
%           Printed: peak depth error, -6 dB lateral width and peak amplitude per target.
% Change  : c1, layer thickness d1, target depths, Fnum.
% Expect  : phase-shift peaks within a few um of the truth, lateral width 0.555 mm at
%           every depth. Stolt at 1540 m/s: every target 0.50 mm too deep (theory
%           (c2/c1 - 1) d1 = 0.497 mm) and ~15 % wider. Stolt at cbar = 1510 m/s: depth exact
%           at 25 mm only (+0.20 mm at 15 mm, -0.20 mm at 35 mm).
% Model   : Born, point elements, no attenuation, no transmission losses or multiples at
%           the interface, IQ at 20 MHz, pitch 0.30 mm, fan filter F# 1.5.
%           Runtime ~3 s in Octave.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);
tStart = tic;

%% Parameters
c2 = 1.540; c1 = 1.450; d1 = 8;               % mm/us, mm
f0 = 5; FBW = 0.6; fs = 20; lam = c2/f0;
Nel = 128; pitch = 0.30; xe = ((0:Nel-1) - (Nel-1)/2)*pitch;
Fnum = 1.5; alphaMax = atan(1/(2*Fnum));
ptsX = [0 0 0 6].';  ptsZ = [15 25 35 25].';  nt = 1200;
cz = @(z) c1*(z < d1) + c2*(z >= d1);          % speed profile
sig_f = FBW*f0/(2*sqrt(2*log(2)));  sig_t = 1/(2*pi*sig_f);
kt = (-ceil(4*sig_t*fs):ceil(4*sig_t*fs)).';  envK = exp(-(kt/fs).^2/(2*sig_t^2));

%% Channel data: vertical transmit, Fermat receive path through the interface
Ttx = d1/c1 + (ptsZ - d1)/c2;                  % us, exact for a normal plane wave
xa = repmat(xe, numel(ptsX), 1);  xb = repmat(ptsX, 1, Nel);   % bracket for the crossing point
lo = min(xa, xb);  hi = max(xa, xb);  Zt = repmat(ptsZ, 1, Nel);  Xt = xb;  Xe = xa;
for it = 1:60                                  % bisection on dT/dxr (monotonic in xr)
    xr = (lo + hi)/2;
    g  = (xr - Xe)./(c1*sqrt((xr - Xe).^2 + d1^2)) - (Xt - xr)./(c2*sqrt((Xt - xr).^2 + (Zt - d1).^2));
    hi(g > 0) = xr(g > 0);  lo(g <= 0) = xr(g <= 0);
end
Trx = sqrt((xr - Xe).^2 + d1^2)/c1 + sqrt((Xt - xr).^2 + (Zt - d1).^2)/c2;
tau = Ttx + Trx;                               % two-way delay (us), targets x elements
q = tau*fs;  i0 = floor(q);  fr = q - i0;  w = exp(-1i*2*pi*f0*tau);
li = i0 + 1 + nt*(0:Nel-1);
Sd = accumarray(li(:), w(:).*(1-fr(:)), [nt*Nel 1]) + accumarray(li(:)+1, w(:).*fr(:), [nt*Nel 1]);
S  = conv2(reshape(Sd, nt, Nel), envK, 'same');
% Snell check: the ray to the edge element of the 15 mm target
[~, k] = max(abs(xe));
s1 = (xr(1,k) - xe(k))/sqrt((xr(1,k) - xe(k))^2 + d1^2);  s2 = (ptsX(1) - xr(1,k))/sqrt((ptsX(1) - xr(1,k))^2 + (ptsZ(1) - d1)^2);
fprintf('Snell check (15 mm target, edge element): sin(a1)/c1 = %.5f, sin(a2)/c2 = %.5f (1/(mm/us))\n', s1/c1, s2/c2);

%% Common spectra and grids
nfft = 2048;  Df = fs/nfft;  fq = [0:nfft/2-1, -nfft/2:-1]*Df + f0;
bq = find(fq >= 0.8 & fq <= 9.2);  [fB, o] = sort(fq(bq));  bq = bq(o);  fB = fB(:);  nB = numel(fB);
nxp = 256;  kx = [0:nxp/2-1, -nxp/2:-1]/(nxp*pitch);  xs = xe(1) + (0:nxp-1)*pitch;
D = fft(S, nfft, 1);  D = fft(D(bq, :), nxp, 2);
D = D .* (abs(kx) <= fB*sin(alphaMax)/c2);     % receive fan, F# 1.5 in tissue
dz = 0.0385;  zg = (0:round(40/dz)).'*dz;  Nz = numel(zg);

%% Phase-shift migration through the layered profile
t1 = tic;  psi = D;  imPS = zeros(Nz, nxp);
stp = cell(1,2);                               % one depth step in each medium
for m = 1:2
    cl = c1*(m == 1) + c2*(m == 2);
    kzr = sqrt(max((fB/cl).^2 - kx.^2, 0));  prop = (fB/cl).^2 >= kx.^2;
    stp{m} = exp(1i*2*pi*(kzr + fB/cl)*dz) .* prop;   % receive kz + transmit f/c
end
for k = 1:Nz
    imPS(k, :) = sum(psi, 1);                  % imaging condition: sum over f
    psi = psi .* stp{1 + (cz(zg(k) + dz/2) == c2)};
end
imPS = ifft(imPS, [], 2);  tPS = toc(t1);

%% Stolt (exact 0-degree map) with a single speed
cbar = 25/(d1/c1 + (25 - d1)/c2);              % time-average speed to 25 mm
ntq = 4096;  Dq = fs/ntq;  fq2 = [0:ntq/2-1, -ntq/2:-1]*Dq + f0;
b2 = find(fq2 >= 0.8 & fq2 <= 9.2);  [fB2, o] = sort(fq2(b2));  b2 = b2(o);  fB2 = fB2(:);  n2 = numel(fB2);
tc = nt/(2*fs);
D2 = fft(S, ntq, 1);  D2 = fft(D2(b2, :) .* exp(1i*2*pi*fB2*tc), nxp, 2);
NZ = 2048;  Kz = (0:NZ-1).'/(NZ*dz);  [KX, KZ] = meshgrid(kx, Kz);
imS = cell(1,2);  cs = [c2 cbar];
for m = 1:2
    c = cs(m);  Dm = D2 .* (abs(kx) <= fB2*sin(alphaMax)/c);
    F = c*(KX.^2 + KZ.^2)./(2*KZ + eps);  J = abs(c*(KZ.^2 - KX.^2)./(2*KZ.^2 + eps));
    ok = KZ > 0 & (KZ - F/c) >= 0;  qq = (F - fB2(1))/Dq + 1;  k0 = floor(qq);  ww = qq - k0;
    ok = ok & k0 >= 1 & k0 < n2;  lk = k0(ok) + n2*(ceil(find(ok)/NZ) - 1);
    I = zeros(NZ, nxp);  I(ok) = ((1 - ww(ok)).*Dm(lk) + ww(ok).*Dm(lk+1)) .* J(ok) .* exp(-1i*2*pi*F(ok)*tc);
    im = ifft2(I)/(dz*Df);  imS{m} = im(1:Nz, :);   % same scale as the phase-shift sum over f bins
end

%% Metrics
names = {'phase-shift', 'Stolt 1540', sprintf('Stolt %.0f', cbar*1e3)};
ims = {imPS, imS{1}, imS{2}};
fprintf('time-average speed to 25 mm: %.1f m/s; expected depth error at 1540 m/s for a target at z:\n', cbar*1e3);
fprintf('   dz = (c2/c1 - 1) d1 = %.3f mm (one-way time preserved), same for every target below the layer\n\n', (c2/c1 - 1)*d1);
fprintf('%-12s %8s %10s %9s %9s\n', 'method', 'target', 'dz peak', 'lat-6dB', 'peak');
fprintf('%-12s %8s %10s %9s %9s\n', '', '(x,z)', '(um)', '(mm)', '(dB)');
xq = -1.5:0.005:1.5;  refpk = zeros(1, numel(ptsX));
for m = 1:3
    E = abs(ims{m});  if m == 1, nrm = max(E(:)); end;  E = E/nrm;
    for it = 1:numel(ptsX)
        sx = abs(xs - ptsX(it)) <= 3;  sz = abs(zg - ptsZ(it)) <= 2;
        sub = E(sz, sx);  [pk, k] = max(sub(:));  [kz, kx2] = ind2sub(size(sub), k);
        zz = zg(sz);  xx = xs(sx);
        zfine = zz(kz) + (-0.2:0.001:0.2);  Lc = interp1(zz, sub(:, kx2), zfine, 'spline');  [pk, b] = max(Lc);
        rowv = interp1(zz, sub, zfine(b), 'spline');  Lr = interp1(xx, rowv, xx(kx2) + xq, 'spline');
        wl = sum(Lr >= max(Lr)/2)*0.005;
        fprintf('%-12s %8s %10.0f %9.3f %9.1f\n', names{m}, sprintf('%g,%g', ptsX(it), ptsZ(it)), (zfine(b) - ptsZ(it))*1e3, wl, 20*log10(pk));
    end
end
fprintf('\nphase-shift loop: %.2f s for %d depth steps (cost ~ Nz Nf Nkx); a Stolt pass: one interpolation\n', tPS, Nz);
fprintf('total run time %.1f s\n', toc(tStart));

%% Figures
figure('Position', [50 50 1100 420]);
for m = 1:3
    subplot(1,3,m); E = abs(ims{m});
    imagesc(xs, zg, max(20*log10(E/max(E(:))), -50), [-50 0]); axis image; axis([-12 12 5 40]); colormap(gray); hold on;
    plot([-12 12], [d1 d1], 'y--'); plot(ptsX, ptsZ, 'r+');
    title(sprintf('%s (dB); red + = truth', names{m})); xlabel('x (mm)'); ylabel('z (mm)');
end

What to look for: the Snell check prints equal \(\sin\alpha/c\) on both sides; phase shift within 1 µm of every target at 0.555 mm width; Stolt at 1540 m/s +0.499 mm (theory 0.497 mm) and 0.635–0.645 mm wide; Stolt at 1510 m/s exact at 25 mm, ±0.20 mm at 15 and 35 mm.

24.12Self-check

Q1At 7.5 MHz, where is the cone edge in \(k_x\), and which arrivals alias with a 0.30 mm pitch?

\(f/c = 7.5/1.54 = 4.87\) cycles/mm. Nyquist is \(1/(2 \times 0.30) = 1.67\) /mm, so arrivals with \(\sin\alpha > 1.67/4.87 = 0.342\), i.e. α > 20.0°, alias.

Q2Zero-offset data of a point at 20 mm are migrated with \(v = c\). Where does it appear, and is it focused?

At 40 mm: \(k_z \approx f/c\) instead of \(2f/c\) doubles depth. It is not focused, because the hyperbola curvature corresponds to the wrong speed; the lab shows 21–24 dB peak loss for the same error at 10 mm.

Q3Compute Garcia's \(\hat{v}/c\), β and γ at θ = 15°.

\(\cos 15° = 0.9659\), \(\sin^2 15° = 0.0670\), so \(1 + \cos\theta + \sin^2\theta = 2.0329\). \(\hat{v}/c = 1/\sqrt{2.0329} = 0.7014\); \(\beta = 1.9659^{1.5}/2.0329\) = \(2.7564/2.0329 = 1.356\); \(\gamma = 0.2588/(2 - 0.9659) = 0.250\).

Q4A 60 µs record is migrated with linear interpolation. What bin spacing Δf keeps the smear at −30 dB, with and without time-centring? FFT length at 20 MHz?

Need \(\varphi = 2\pi\,\Delta f\,\lvert t - t_c\rvert \le 0.92\), i.e. \(\Delta f\,\lvert t - t_c\rvert \le 0.146\). Uncentred, \(\lvert t\rvert \le 60\) µs: \(\Delta f \le 2.44\) kHz, 8192 points. Centred, \(\lvert t - t_c\rvert \le 30\) µs: \(\Delta f \le 4.88\) kHz, 4096 points.

Q5Estimate the narrowband lateral −6 dB width at 5 MHz for a continuous plane-wave span of ±10° with receive F# 1.5.

Lateral k-space is the convolution of the transmit span (half-width \((f/c)\sin 10° = 0.564\) /mm) and the receive fan (\((f/c)\sin\alpha = 1.027\) /mm), so the PSF is \(\mathrm{sinc}(2 \times 0.564\, x)\cdot\mathrm{sinc}(2 \times 1.027\, x)\). Its −6 dB width is 0.53 mm, against 0.59 mm receive-only. A discrete 7-angle set at ±9° simulated 0.505 mm: a discrete transmit factor has a slightly narrower central lobe.

Q6A 6 mm fat layer (1450 m/s) lies over tissue (1540 m/s). How far off is a single-speed (1540 m/s) Stolt image below it?

\(\Delta z = (c_2/c_1 - 1)\,d_1\) = \((1.0621 - 1) \times 6 = 0.37\) mm too deep, the same for every target below the layer. Phase-shift migration with the layer speed removes it.

Q7Why does compounding N angles in f-k need only one inverse 2D FFT, and what does it not fix?

The inverse FFT is linear: \(\sum_\theta \mathrm{IFFT}(R_\theta) = \mathrm{IFFT}(\sum_\theta R_\theta)\), so remapped spectra are summed on the common grid first. Each angle still needs its own forward FFTs, steering compensation and interpolation, and motion between angles still adds phase errors.

24.13Figures from the literature

Published figures with real or carefully simulated data that show what this module computes. Links only (no reproduced images); figure numbers checked against the cited version.

Fourier-domain plane-wave imaging against delay-and-sum on a wire phantom

J. Cheng and J.-y. Lu, "Extended high-frame rate imaging method with limited-diffraction beams," IEEE TUFFC 53(5), 880–899, 2006. Fig. 8 (steered plane waves; Fig. 7 is the limited-diffraction-beam version). doi:10.1109/TUFFC.2006.1632680

Wires in water reconstructed in the Fourier domain from 1, 11, 19 and 91 steered plane waves, beside delay-and-sum with a fixed 60 mm transmit focus and with a synthesized dynamic transmit focus. Compare the Fourier reconstructions as transmissions are added with the two DAS panels: more angles means more crescents of Figure 24.B filling the coverage. This is the k-space map of (24.14) on experimental data.

The same method in vivo: kidney and heart

Same source. Figs. 13 and 14. doi:10.1109/TUFFC.2006.1632680

Fig. 13: a human right kidney from 91 steered plane waves against delay-and-sum with a 70 mm fixed focus and 88 transmissions. Fig. 14: an ECG-triggered heart with 11, 19 and 91 transmissions. Compare the image quality at equal transmit count: the Fourier method needs no per-pixel delay loop, and its quality is set by the angle coverage, as in Demo 24.4.

Stolt and Lu f-k against time-domain plane-wave imaging, with timings

L. Merabet, S. Robert and C. Prada, "2-D and 3-D reconstruction algorithms in the Fourier domain for plane-wave imaging in nondestructive testing," IEEE TUFFC 66(4), 772–788, 2019. Figs. 6–8 (2-D images, echodynamic curves, close-ups), Fig. 9 (computation time against pixel count), Figs. 16–17 (experimental 3-D). doi:10.1109/TUFFC.2019.2895995

Side-drilled hole pairs 2 mm apart in steel imaged by Stolt, Lu and time-domain plane-wave imaging (PWI): the three agree in image quality under the aperture (outside it Stolt blurs and Lu is damped at high angles, Figs. 7–8), which is the "same physics, different bookkeeping" of Section 24.5 in a different medium. Fig. 9 is the cost argument: the Fourier methods scale as N log N. The text reports, for 1000 × 1000 pixels in MATLAB, Lu 2.05 s, Stolt 3.8 s and PWI 41.4 s, and for a 3-D volume of 100 × 100 × 490 voxels with 25 plane waves, Lu 129.6 s, Stolt 136.8 s, PWI 1 h 18 min (plotted below). In the experimental 3-D case PWI took 38 times longer than Lu, against a predicted factor of 36.

Approx. modelRe-plotted from the timings published in the text of Merabet et al. (2019): their MATLAB code and hardware, not re-measured here.

Matrix imaging: what lies beyond the confocal f-k image

W. Lambert, L. A. Cobus, T. Frappart, M. Fink and A. Aubry, "Distortion matrix approach for ultrasound imaging of random scattering media," PNAS 117(26), 14645–14656, 2020; preprint arXiv:1912.05417. Fig. 2 (numbering checked in the arXiv version). doi:10.1073/pnas.1921533117 · arXiv

A reflection matrix acquired through a plexiglass layer over a tissue-mimicking phantom: the original confocal image, the image after multiple reflections are removed, and the image after matrix aberration correction, each with a map of the Strehl ratio. The confocal image is what Section 24.4 forms; the correction uses the off-confocal terms that f-k throws away (Section 24.6), which is why matrix methods are the natural next step for aberration correction.

24.14Go deeper

R. H. Stolt, "Migration by Fourier transform," Geophysics 43(1), 23–48, 1978. doi:10.1190/1.1440826. The change of variables \(f \to k_z\) and the interpolation it requires.
J. Gazdag, "Wave equation migration with the phase-shift method," Geophysics 43(7), 1342–1351, 1978. doi:10.1190/1.1440899. Depth stepping in the angular spectrum; the layered extension.
D. Loewenthal, L. Lu, R. Roberson and J. Sherwood, "The wave equation applied to migration," Geophysical Prospecting 24, 1976. The exploding-reflector model.
J. F. Claerbout, Imaging the Earth's Interior, Blackwell, 1985. Migration from first principles.
P. L. Stoffa, J. T. Fokkema, R. M. de Luna Freire and W. P. Kessinger, "Split-step Fourier migration," Geophysics 55(4), 1990. Mild lateral velocity variation in the Fourier domain.
D. Garcia, L. Le Tarnec, S. Muth, E. Montagnon, J. Porée and G. Cloutier, "Stolt's f-k migration for plane wave ultrasound imaging," IEEE TUFFC 60(9), 1853–1867, 2013. doi:10.1109/TUFFC.2013.2771. The ERM fit derived in Section 24.3.
J. Cheng and J.-y. Lu, "Extended high-frame rate imaging method with limited-diffraction beams," IEEE TUFFC 53(5), 880–899, 2006 (doi); J.-y. Lu, "2D and 3D high frame rate imaging with limited diffraction beams," IEEE TUFFC 44(4), 1997. Fourier-domain plane-wave and limited-diffraction-beam imaging.
L. Merabet, S. Robert and C. Prada, IEEE TUFFC 66(4), 772–788, 2019 (doi). Stolt and Lu in 2-D and 3-D, with complexity analysis.
A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, 1988, chapter 6. The Fourier diffraction theorem behind the Ewald-circle map.
G. Montaldo, M. Tanter, J. Bercoff, N. Benech and M. Fink, "Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography," IEEE TUFFC 56(3), 489–506, 2009. doi:10.1109/TUFFC.2009.1067.
N. Bottenus, "Recovery of the complete data set from focused transmit beams," IEEE TUFFC 65(1), 30–38, 2018 (REFoCUS). doi:10.1109/TUFFC.2017.2773495.
W. Lambert, L. A. Cobus, T. Frappart, M. Fink and A. Aubry, "Distortion matrix approach for ultrasound imaging of random scattering media," PNAS 117(26), 14645–14656, 2020. doi:10.1073/pnas.1921533117.
T. Chernyakova and Y. C. Eldar, "Fourier-domain beamforming: the path to compressed ultrasound imaging," IEEE TUFFC 61(8), 2014; preprint arXiv:1307.6345. The other meaning of the term.
J. A. Jensen, S. I. Nikolov, K. L. Gammelmark and M. H. Pedersen, "Synthetic aperture ultrasound imaging," Ultrasonics 44, 2006. Time-domain synthetic aperture for comparison.

← 23 Beams, k-space, compounding 25 Adaptive & coherence BF →

24.1The angular spectrum: channel data as a sum of plane waves

Channel data and its kx–f spectrum: cone, Nyquist and folding

24.2Exploding reflectors and Stolt's change of variables

Stolt mapping, step by step (zero-offset data, exploding reflector)

24.3Plane-wave transmits: steering compensation and Garcia's ERM fit

The ERM fit: true travel time against the fitted hyperbola

24.4The exact map: transmit plus receive wavevectors

k-space coverage of a plane-wave set and the PSF it implies

24.5f-k against DAS: same physics, different bookkeeping

f-k against DAS on the same channel data: images, difference, timing

24.6Beyond one speed of sound: layers, sound-speed sweeps and the reflection matrix

24.7Visual summary

24.8How experts think about it

24.9Cheat sheet

24.10Interview and qualifying-exam questions

24.11MATLAB lab

24.12Self-check

24.13Figures from the literature

Fourier-domain plane-wave imaging against delay-and-sum on a wire phantom

The same method in vivo: kidney and heart

Stolt and Lu f-k against time-domain plane-wave imaging, with timings

Matrix imaging: what lies beyond the confocal f-k image

24.14Go deeper

Channel data and its k_x–f spectrum: cone, Nyquist and folding