Simulation and characterization tools: spatial impulse responses, full-wave codes, and the measurement bench
Every claim about a probe or a beamformer ends in a simulation or a measurement. This module rebuilds the linear theory under Field II and FOCUS (the spatial impulse response), shows what full-wave codes such as k-Wave, FDTD and explicit finite elements add and what they cost, and then turns to the bench: pulse-echo tests, impedance, hydrophone scans, acoustic output, array QA, and how measurements are turned back into a model.
You will be able to
- Compute and read the spatial impulse response of a rectangular element on and off axis, and predict the field it radiates
- Choose between a linear field simulator and a full-wave code for a given question, and estimate the grid, points per wavelength and memory the full-wave run needs
- Recognize numerical dispersion and other simulator errors, and validate a code against analytic benchmarks
- Extract centre frequency, −6 dB bandwidth, pulse length and loop sensitivity from a pulse-echo waveform without the usual definition errors
- Correct hydrophone data for spatial averaging and bandwidth, and connect scans to MI, TI and the IEC and FDA documents
- Quantify dead elements and crosstalk, and turn bench data into KLM and FE model parameters
Assumes
- KLM and impedance of a plate (Module 16, Piezoelectrics & KLM); matching, lens and FE design (Module 17, Matching, backing, lens, FEA)
- Element directivity and array factors (Module 19, Arrays: 1D to 2D, RCA, CMUT); front-end chain (Module 20, Front-end electronics)
- Aperture ↔ beam Fourier pair (Module 04, Sampling, Fourier, LTI); analytic signal and envelope (Module 08, Analytic signal & IQ)
- Wave equation, diffraction and Fubini nonlinearity (Module 15, Wave physics & harmonics); MI, TI and derating are developed in Module 22 (Acoustic output & safety)
21.1Two families of simulators
Ultrasound simulators split along one line: do they assume the medium is linear and homogeneous, or do they solve the wave equation over a grid of the actual medium?
Linear field simulators (Field II, FOCUS) assume a linear, homogeneous medium with a single sound speed. The field of any aperture is then a superposition of free-space Green's functions over its surface, which collapses into one function of position and time: the spatial impulse response \(h(\mathbf{r}, t)\). Scattering is added as a sum of point scatterers in the Born approximation, so an image with 10⁵ to 10⁶ scatterers is routine. They cannot represent a sound-speed change, refraction, reverberation between structures, or harmonic generation.
Full-wave simulators (k-Wave, FDTD codes, explicit finite elements such as PZFlex) discretize the wave equations over the whole domain. They handle heterogeneity, multiple scattering, power-law absorption and nonlinearity, and finite elements also model the transducer structure itself. The price is a grid that must resolve the shortest wavelength everywhere, which in 3D grows as \(f^4\) in run time.
| Tool | Method | Models well | Cannot model or costly | Typical use |
|---|---|---|---|---|
| Field II (Jensen 1996; Jensen & Svendsen 1992) | Spatial impulse response; aperture split into small rectangles or triangles, each treated in its own far field; pulse-echo by convolving transmit and receive responses | Arbitrary apertures, apodization, focusing, pulsed fields, point-scatterer phantoms, flow | Heterogeneous sound speed, multiple scattering, nonlinearity; attenuation only as a linearized frequency-dependent term | Beam design, B-mode and Doppler image simulation, beamformer testing |
| FOCUS (McGough group, Michigan State) | Fast nearfield method (removes the singularity of the Rayleigh integral near the face), time-space decomposition for transients, angular spectrum for propagation | Accurate near fields of pistons and arrays, CW and transient, HIFU arrays | Same linear, homogeneous limits as Field II | Reference near-field calculations, therapy array design |
| k-Wave (Treeby & Cox 2010) | k-space pseudospectral time domain on a regular grid, PML boundaries | Heterogeneous c and ρ, power-law absorption, nonlinear propagation, an elastic variant; GPU/C++ versions | Memory and time grow with \((\text{domain}/\lambda)^3\); staircased geometry; the transducer is a source distribution, not a structure | Aberration, skull and fat layers, sound-speed estimation tests, nonlinear and therapy fields |
| FDTD (generic) | Finite differences on a staggered grid, leapfrog in time | Same physics as k-Wave; simple to write and extend | Needs many more points per wavelength (numerical dispersion) | Teaching, custom physics, 2D studies |
| Explicit FE: PZFlex (now in the OnScale / Ansys family) | Time-domain finite elements with piezoelectric elements and an electrical port | Kerfs, layers, lens, lateral modes, crosstalk, impedance and pulse-echo of a real stack | Large tissue domains are expensive; usually coupled to an extrapolation or a linear code for the far field | Transducer and array design (Module 17) |
| Implicit or harmonic FE (e.g. COMSOL) | Frequency-domain solves, eigenmodes | Impedance spectra, mode shapes of small structures, CMUT cells | Broadband transients need many frequency solves | Element and cell design, mode identification |
What actually runs on this page. Every simulated curve below comes from code written for this course and executed in your browser or in the Octave scripts. The commercial and research packages are described from their papers and manuals; a name such as "Field-II-style" in a demo title only says which idea the demo reimplements.
| Solver or tool | Status here | Implementation | Where |
|---|---|---|---|
| Exact spatial impulse response of a rectangle (arc method) | Run | own JS; own Octave | Demo 21.1, m21_spatial_impulse_response.m |
| Field-II-style sum of far-field sub-elements | Run (reimplementation of the idea, not Field II) | own JS; own Octave | Demo 21.2, same script |
| 2D staggered-grid FDTD, FD2 and FD4 in space, absorbing sponge | Run | own JS, on the button | Demo 21.4 |
| 2D linear model (sum of 2D Green's functions) | Run | own JS | Demo 21.4 |
| 1D FD2, FD4 and k-space pseudospectral time stepping | Run | own Octave | m21_numerical_dispersion.m |
| Discrete dispersion relations, Fubini series, disc averaging, BVD impedance, Rayleigh integral of a focused disc | Evaluated | closed forms and quadrature, JS and Octave | Demos 21.3 and 21.6, cards, rule figures |
| Field II, FOCUS, k-Wave, PZFlex/OnScale, COMSOL | Described only | not run here | tables in 21.1 and 21.5; published comparisons in 21.15 |
ULM with linear microbubble scatterers and a homogeneous medium: a linear simulator is enough and fast (millions of point scatterers over thousands of frames). Sound-speed estimation and aberration correction: the effect being studied is the heterogeneity itself, so the forward model must be full-wave (k-Wave or FDTD) or at least include the layer explicitly; a linear simulator with a "wrong" global c only tests the global-speed problem. Nonlinear bubble response or harmonic imaging: a nonlinear full-wave code or a bubble model, never Field II alone.
21.2The spatial impulse response
A planar aperture S in a rigid baffle vibrates with uniform normal velocity \(v(t)\). The Rayleigh integral gives the velocity potential at \(\mathbf{r}\) as a sum of delayed, spherically spreading copies of \(v\) from every surface point. Because every copy is the same \(v(t)\), the geometry factors out into one convolution:
- \(\varphi\)
- velocity potential, m²/s
- \(\mathbf{r},\ \mathbf{r}'\)
- field point and source point on the aperture, m
- \(S\)
- radiating aperture (in a rigid baffle)
- \(v(t)\)
- uniform normal surface velocity, m/s
- \(c\)
- sound speed, m/s
- \(h(\mathbf{r}, t)\)
- spatial impulse response, m/s
- \(\rho_0\)
- ambient density, kg/m³
- \(p\)
- pressure, Pa
- \(*\)
- convolution in time
- \(t\)
- time, s
- \(\delta\)
- Dirac delta function, 1/s
- \(dS\)
- surface element of the aperture, m²
- \(\partial\)
- partial derivative
The geometric reading. Let \(z\) be the distance from the aperture plane and \(\sigma\) the radius in that plane around the projection of \(\mathbf{r}\). Then \(R^2 = \sigma^2 + z^2\), so \(R\, dR = \sigma\, d\sigma\), and \(dS = \sigma\, d\sigma\, d\theta = R\, dR\, d\theta\). The delta picks \(R = ct\), and \(\delta(t - R/c)\, dR\) integrates to \(c\):
- \(z\)
- distance of the field point from the aperture plane, m
- \(\sigma\)
- radius, in the aperture plane, of the circle around the projection of \(\mathbf{r}\), m
- \(\Theta\)
- total arc angle of that circle inside the aperture, rad
- \(t\)
- time after the pulse, s
- \(h\)
- spatial impulse response, as in (21.1)
Disc, on axis. The circle stays inside the disc of radius \(a\) from \(t = z/c\) to \(t = \sqrt{z^2 + a^2}/c\), so \(h\) is a rectangle of height \(c\), and the pressure is two copies of the source waveform: \(p = \rho_0 c\, [v(t - z/c) - v(t - \sqrt{z^2 + a^2}/c)]\). The first is the plane wave from the face, the second an inverted edge wave. Near the face they are separated in time; far away they overlap and the difference becomes a derivative, which is the far-field point-source limit.
Rectangle. With half-widths \(a \le b\) on axis, \(\Theta = 2\pi\) until \(\sigma\) reaches \(a\), then each pair of edges cuts off \(4 \arccos(\cdot/\sigma)\), and \(h\) ends at the corners:
- \(a,\ b\)
- half-widths of the rectangle, \(a \le b\), m
- \([\,\cdot\,]\)
- indicator: 1 when the condition holds, else 0
- \(\sigma,\ \Theta\)
- as in (21.2)
- \(\arccos\)
- inverse cosine, rad
- \(y(t)\)
- received echo voltage, V
- \(E_m(t)\)
- electromechanical impulse response of the transducer
- \(v\)
- surface velocity waveform, m/s
- \(h_t,\ h_r\)
- transmit and receive spatial impulse responses, m/s
- \(\propto\)
- proportional to
- \(t,\ \mathbf{r}\)
- time (s) and field point (m)
- \(*\)
- convolution in time
- \(\partial^3/\partial t^3\)
- third time derivative of the velocity: one from radiation, two from Rayleigh scattering by a small inhomogeneity
Spatial impulse response of a rectangular element, on and off axis
SimulationOwn JS: \(\Theta(t)\) found by testing points on each arrival circle against the rectangle, then \(p = \rho_0 v' * h\) by FFT. The same method runs in m21_spatial_impulse_response.m (Octave). Not Field II output.For a receive beamformer you ask "when does the echo from this pixel reach each element?" The SIR asks the same question for every patch of one element and histograms the answers. A field point that sees all patches at nearly the same delay (far field, on axis) gets a short, tall response; a point with a wide spread of delays (near the face, or far off axis) gets a long, low one that low-pass filters the pulse.
Exact SIRs jump (the on-axis disc response is a rectangle). Sampling them directly at a typical simulation rate puts energy errors on every edge. Field II therefore tracks the exact break times and works with integrated responses; its user guide recommends sampling around 100 MHz with these approximations rather than GHz rates with exact expressions. If you write your own SIR code, integrate over each sample interval instead of point-sampling.
21.3From the impulse response to the field: Field II and FOCUS
The continuous-wave field is one Fourier coefficient of the same function: \(P(\mathbf{r}, f) = j2\pi f \rho_0 V(f) H(\mathbf{r}, f)\). Everything a linear code computes (CW beam plots, pulsed peak pressures, pulse-intensity integrals, pulse-echo signals) is a convolution or a Fourier transform of \(h\).
Field II's approximation. Exact \(h\) is known in closed form only for simple shapes, so Field II splits the aperture into small rectangles and uses a far-field response for each: a short pulse of area \(\Delta A/(2\pi R_i)\) at \(t = R_i/c\). The sum converges to the exact \(h\) when each sub-element is in its own far field, roughly \(R \gg d^2/\lambda\) for sub-element size \(d\). Demo 21.2 compares the subdivided field with the exact on-axis response so you can see where that breaks.
FOCUS's approach. The fast nearfield method rewrites the Rayleigh integral over a rectangle or disc as a one-dimensional integral along its edges with a smooth integrand, which removes the \(1/R\) singularity that makes direct sums inaccurate close to the face. The time-space decomposition gives transient fields at lower sampling rates, and the angular spectrum approach propagates a field from one plane to the next. Same physics as Field II, with better accuracy near the transducer per unit of computation.
- \(v_0\)
- surface velocity amplitude, m/s
- \(a\)
- disc radius, or half the side of the square, m
- \(\lambda\)
- wavelength, m
- \(z\)
- distance along the axis, m
- \(p(z)\)
- on-axis pressure amplitude, Pa
- \(\rho_0,\ c\)
- density (kg/m³) and sound speed (m/s)
- \(\lvert\cdot\rvert\)
- magnitude
The field of a rectangular element: subdivided sum versus the exact response
SimulationOwn JS reimplementation of the Field II idea (point responses \(\Delta A/(2\pi R)\) of \(n \times n\) sub-elements) against the closed-form on-axis SIR; "pulsed" is the rms over a 60 % Gaussian band. Field II itself was not run.21.4Full-wave simulators: FDTD, k-space pseudospectral and explicit finite elements
All three solve the same first-order system, here in its linear lossless form with a mass source \(q\):
- \(\mathbf{v}\)
- particle velocity, m/s
- \(p\)
- acoustic pressure, Pa
- \(q\)
- mass source (volume injection rate per unit volume), 1/s
- \(\rho_0,\ c\)
- density (kg/m³) and sound speed (m/s), both may vary with position
- \(t\)
- time, s
- \(\nabla,\ \nabla\cdot\)
- gradient and divergence, 1/m
FDTD. Pressure and velocity live on staggered grids and leapfrog in time. Every scheme of this type has a dispersion relation linking the numerical frequency \(\omega\) to the wavenumber \(k\). With Courant number \(\mathrm{CFL} = c\Delta t/\Delta x\):
- \(\omega\)
- numerical angular frequency, rad/s
- \(k\)
- wavenumber, rad/m
- \(\Delta x,\ \Delta t\)
- grid spacing (m) and time step (s)
- \(\mathrm{CFL}\)
- Courant number \(c\Delta t/\Delta x\), dimensionless
- \(S\)
- the scheme's discrete spatial-derivative operator, dimensionless
- \(c_{\mathrm{ref}}\)
- reference sound speed of the k-space operator, m/s
- \(j\)
- imaginary unit
- \(\operatorname{sinc}\)
- \(\sin(\pi u)/(\pi u)\)
- \(x\)
- position on the grid, m
k-space pseudospectral (k-Wave). Spatial derivatives are taken with FFTs, which are exact for band-limited fields, so the grid can approach the Nyquist limit of two points per wavelength (a few more in practice, to represent heterogeneities and sources). The temporal error of leapfrog is then cancelled by the k-space correction (Tabei, Mast & Waag, JASA 2002): in a homogeneous medium at \(c_{\mathrm{ref}}\) the time stepping is exact for any \(\Delta t\), and in heterogeneous media the error is small when \(c\) stays near \(c_{\mathrm{ref}}\). Perfectly matched layers absorb waves leaving the grid.
Explicit FE (PZFlex class). The same leapfrog idea on a mesh of elements, with piezoelectric coupling and an electrical port, so a kerfed, layered, lensed stack is modelled as it is built (Module 17). Mesh density is chosen by refinement studies; the step is limited by the smallest element.
Numerical dispersion: why FD2 needs many points per wavelength
Exact formulaEach scheme's discrete dispersion relation; the propagated pulses apply that relation to a Gaussian pulse in the Fourier domain on a periodic grid (equivalent to time stepping without boundaries).Linear versus full-wave on the same problem. Demo 21.4 runs a small 2D full-wave simulation written for this page (staggered FDTD with an absorbing sponge, 0.1 mm grid, 2 MHz, 32 elements at 0.3 mm pitch focused at 16 mm) and a 2D linear model (sum of 2D Green's functions from the same source cells). In a homogeneous medium they must agree; with a fat layer only the full-wave code sees the change.
Linear (Green's-function) model versus a 2D FDTD full-wave run, both own code
Simulation2D staggered-grid FDTD (FD4 or FD2 in space, leapfrog in time, absorbing sponge), 0.1 mm grid, run in JS when you press the button; the linear model sums 2D Green's functions from the same source cells. It plays the role that k-Wave and Field II play in research; neither was run.A focus that moves or widens in a simulation is exactly the signature of a sound-speed error or an aberration, so grid dispersion can masquerade as the effect you are studying. Before trusting any aberration or sound-speed result, rerun the homogeneous case on the same grid and confirm it matches a linear reference; then refine the grid until the quantity you report stops changing.
21.5What each tool can model, and how to validate it
| Phenomenon | Field II / FOCUS | k-Wave / FDTD | Explicit FE (PZFlex class) |
|---|---|---|---|
| Diffraction, focusing, apodization, pulsed fields | Yes (fast) | Yes | Yes (near field; far field by extrapolation) |
| Frequency-dependent attenuation | Approximate (linearized, path-averaged) | Yes (power law) | Yes (material damping models) |
| Sound-speed heterogeneity, refraction, aberration | No | Yes | Yes, at high cost in tissue |
| Multiple scattering, reverberation | No (single scattering) | Yes | Yes |
| Nonlinear harmonics, shocks | No | Yes (k-Wave nonlinear) | Possible in some codes |
| Shear waves, bone, solids | No | Elastic versions | Yes |
| Transducer structure: kerf modes, lens, crosstalk, impedance | No (ideal piston per element) | No (source distribution) | Yes |
| Speckle image from 10⁵–10⁶ scatterers | Yes (routine) | Expensive | Impractical |
Validation ladder. (1) Analytic benchmarks: the on-axis disc formula, the far-field directivity sinc and jinc, the Fubini solution for weak plane-wave nonlinearity. (2) Convergence: halve \(\Delta x\) (and \(\Delta t\)) and confirm the quantity you report changes less than its tolerance; the order of the scheme predicts how fast (\(\Delta x^2\) for FD2; \(\Delta x^4\) for the spatial part of FD4, until the second-order leapfrog time error dominates at fixed CFL). (3) Code to code: a homogeneous-water case in Field II or FOCUS against the full-wave code. (4) Measurement: hydrophone scans of the real probe (21.7), pulse-echo and impedance (21.6), with the hydrophone's own averaging and bandwidth modelled before comparing.
21.6Pulse-echo testing and electrical impedance
Setup. A pulser-receiver excites the transducer in degassed water; the echo from a flat, polished reflector (steel or fused quartz) at normal incidence and at the focal or natural-focus distance is recorded on a scope sampling well above the band (here 250 MS/s). The probe is tilted until the echo is maximal. The gated echo is the transducer's two-way impulse response times the reflector, water and diffraction losses.
What is reported. Peak frequency \(f_{\mathrm{pk}}\); the −6 dB points \(f_l\) and \(f_u\) of the echo amplitude spectrum; centre frequency \(f_c = (f_l + f_u)/2\) (the arithmetic mean, as in ASTM E1065; some datasheets use \(\sqrt{f_l f_u}\)); fractional bandwidth \((f_u - f_l)/f_c\); −6 dB and −20 dB pulse durations from the envelope; and loop sensitivity \(20 \log_{10}(V_{\mathrm{echo,pp}}/V_{\mathrm{tx,pp}})\).
- \(V_{\mathrm{echo,pp}},\ V_{\mathrm{tx,pp}}\)
- peak-to-peak echo and excitation voltages, V
- \(R_{\mathrm{refl}}\)
- amplitude reflection coefficient of the reflector, dimensionless
- \(\alpha_w\)
- attenuation of water, dB/(mm·MHz²)
- \(z\)
- stand-off distance to the reflector, mm
- \(f_c\)
- centre frequency, MHz
- \(R_{\mathrm{steel/water}}\)
- amplitude reflection coefficient of a steel block in water (impedances 45.4 and 1.48 MRayl)
- \(\log_{10}\)
- base-10 logarithm
Pulse-echo test analysis: centre frequency, −6 dB bandwidth, pulse length, loop sensitivity
Approx. modelSynthetic "measured" echo: causal gamma-tone pulse, one lens reverberation and white noise. The analysis (gate, FFT, −6 dB points, envelope, loop sensitivity) is the real procedure.Impedance analyzer. Measured in air and then in water, \(\lvert Z(f) \rvert\) and \(\angle Z(f)\) show the series resonance \(f_r\) (impedance minimum, conductance maximum) and the antiresonance \(f_a\) (impedance maximum). The IEEE Std 176 relation gives the thickness coupling \(k_t^2 = (\pi/2)(f_r/f_a) \tan[(\pi/2)(f_a - f_r)/f_a]\) (Module 16); the capacitance well above resonance (or the low-frequency capacitance times \(1 - k_t^2\)) gives the clamped permittivity \(\varepsilon^S\); the peak sharpness gives the loss. Per element, an open (broken lead) shows almost no capacitance and no resonance, a short shows a near-zero \(\lvert Z \rvert\), and a depoled element keeps its capacitance but loses the resonance.
21.7Hydrophone measurements and acoustic output
Hydrophones. Membrane hydrophones (PVDF film with a small poled spot) have the flattest frequency response and serve as references; needle hydrophones are small and robust but have a less flat response, and at low frequency their effective radius can exceed the geometric one by a large factor (Wear, Fury & Alvarenga, IEEE TUFFC 2025); fibre-optic hydrophones are small and survive therapy pressures. The calibration is the end-of-cable sensitivity \(M_L(f)\) in V/Pa, ideally with phase, so that \(p(t) = \mathcal{F}^{-1}\{V(f)/M_L(f)\}\) can be deconvolved.
The documents, by name. IEC 62127-1, Ultrasonics – Hydrophones – Part 1: Measurement and characterization of medical ultrasonic fields (the 2007 edition was titled "up to 40 MHz"; the 2022 edition removed that limit): how to scan, and the limits on hydrophone size and bandwidth with their corrections. IEC 62127-2, Part 2: Calibration for ultrasonic fields up to 40 MHz, and IEC 62127-3, Part 3: Properties of hydrophones for ultrasonic fields up to 40 MHz. IEC 62359, Ultrasonics – Field characterization – Test methods for the determination of thermal and mechanical indices related to medical diagnostic ultrasonic fields (a newer edition adds "and acoustic intensities" to the title). IEC 60601-2-37, Medical electrical equipment – Part 2-37: Particular requirements for the basic safety and essential performance of ultrasonic medical diagnostic and monitoring equipment. The FDA guidance Marketing Clearance of Diagnostic Ultrasound Systems and Transducers (2019, revised 2023) sets the Track 3 limits: derated \(I_{\mathrm{SPTA.3}}\) ≤ 720 mW/cm², MI ≤ 1.9 (or \(I_{\mathrm{SPPA.3}}\) ≤ 190 W/cm²), TI above 6 to be justified; ophthalmic TI ≤ 1, \(I_{\mathrm{SPTA.3}}\) ≤ 50 mW/cm², MI ≤ 0.23.
- \(p_r\)
- peak rarefactional pressure measured in water, MPa
- \(p_{\mathrm{r.3}}\)
- derated (0.3 dB/cm/MHz) peak rarefactional pressure, MPa
- \(f,\ f_{\mathrm{awf}}\)
- frequency; acoustic working frequency, MHz
- \(z\)
- depth, cm
- \(r\)
- radial distance from the beam axis in the focal plane, mm
- \(F\#\)
- F-number of the circular aperture
- \(b\)
- effective hydrophone radius, mm
- \(\langle p \rangle\)
- pressure averaged over the hydrophone disc
- \(\mathrm{MI}\)
- mechanical index, dimensionless
- \(\lambda\)
- wavelength, mm
- \(\operatorname{jinc}\)
- \(2J_1(\pi u)/(\pi u)\)
- \(p(r),\ p(0)\)
- focal-plane pressure at radius \(r\) and on axis, MPa
Hydrophone spatial averaging at a focus
Approx. modelIdeal circular-aperture focus with a jinc lateral profile, uniform averaging of complex pressure over the hydrophone disc, second harmonic \(\propto p_1^2\) (quasilinear).The acoustic output table in a submission is a chain of corrections: hydrophone calibration (with uncertainty), spatial-averaging and bandwidth corrections, derating at 0.3 dB/cm/MHz, the depth where the derated quantity peaks, and the worst-case operating mode found by a search over settings. Each step can move MI by tens of per cent, which is why the uncertainty budget is part of the report.
A Schlieren system images the refractive-index modulation that the acoustic pressure imposes on water, giving a picture of the whole beam (side lobes, grating lobes, asymmetries) in seconds. It is ideal for quick visual checks and pass/fail screening; quantitative pressure needs calibration against a hydrophone.
21.8Array QA, phantoms, and turning measurements into a model
Element-level tests. Per-element pulse-echo from a flat target gives amplitude, centre frequency and bandwidth maps across the array; per-element impedance identifies opens and shorts. Crosstalk is measured by driving one element and recording the neighbours (electrically, or their surface motion with a laser vibrometer), reported in dB versus neighbour index. Weak or dead elements and crosstalk change the beam in predictable ways, which Demo 21.7 quantifies.
System-level tests. Tissue-mimicking phantoms with a specified sound speed (1540 m/s) and attenuation (commonly 0.5 or 0.7 dB/cm/MHz) carry wire or filament targets for resolution and distance accuracy, and anechoic or hyperechoic cylinders for contrast. IEC 61391-1 (Ultrasonics – Pulse-echo scanners – Part 1: Techniques for calibrating spatial measurement systems and measurement of system point-spread function response) and IEC 61391-2 (Part 2: Measurement of maximum depth of penetration and local dynamic range) standardize these.
- \(\sigma_a\)
- rms relative amplitude error per element, dimensionless
- \(w_n\)
- apodization weight of element n
- \(N,\ D\)
- number of active elements; number of dead elements
- \(\kappa\)
- nearest-neighbour crosstalk amplitude ratio, linear
- \(k\)
- wavenumber, rad/m
- \(p\)
- element pitch, m
- \(\theta\)
- angle from broadside, degrees
- \(\sum\)
- sum over the elements \(n\)
- \(\cos\)
- cosine (argument in rad)
Dead elements, element spread and crosstalk: what they do to the beam
SimulationMonte Carlo array draw in JS: Gaussian amplitude spread in dB, randomly placed dead elements, array factor × rectangular element factor with in-phase nearest-neighbour crosstalk.Closing the loop. Measurements become model parameters in a fixed order, each fit using only the data that constrain it:
| Measurement | Fits in the model | Model |
|---|---|---|
| Impedance of the bare plate in air (\(f_r\), \(f_a\), low-frequency C) | \(k_t\), clamped permittivity \(\varepsilon^S\), mechanical and dielectric loss | KLM / Mason (Module 16) |
| Impedance of the finished element in air and water | Matching-layer and backing impedances and thicknesses, bond lines | KLM with layers (Module 17) |
| Pulse-echo waveform and spectrum | Two-way transfer function, cable and tuning network, receiver load | KLM + electrical network (Module 20) |
| Spurious peaks in impedance, crosstalk, element directivity | Lateral modes, kerf filler, sub-dicing, lens coupling | 2D/3D FE (PZFlex class) |
| Hydrophone field scans | Effective aperture, apodization, lens focal length, element directivity | Field II / FOCUS with fitted parameters |
| Phantom images | System PSF, TGC and dynamic range, sound-speed assumption | Field II image simulation |
21.9Visual 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 in the text (c = 1540 m/s unless stated).
Simulation cards run the arc-method SIR and the Field-II-style subdivision in JS. Model cards: Fubini lossless plane wave; Gaussian transducer bands; two-Gaussian double-humped spectrum; jinc focus; Rayleigh integral of a focused disc for MI; independent element errors.
21.10How experts think about it
An expert first decides which physics the question needs: if the medium can be treated as homogeneous and linear, the spatial impulse response answers everything quickly and exactly; if the answer depends on heterogeneity, nonlinearity or the transducer structure, only a full-wave code will do, and its grid must be proven converged. On the bench, every number is a ratio with a definition behind it (two-way or one-way, −6 dB or −3 dB, arithmetic or geometric centre, derated or not), and every measuring instrument (the hydrophone above all) filters what it measures. Expert habit: model the instrument, then compare.
Linearity and a homogeneous medium make the field a superposition of identical delayed copies of \(v(t)\), so the geometry factors out into \(h(\mathbf{r}, t)\). CW beam plots are \(H\) at one frequency, pulsed maps are \(v' * h\), and pulse-echo is \(h_t * h_r\) with three time derivatives. For a 4 × 6 mm element at z = 15 mm on axis, \(h\) lasts 0.28 µs, shorter than one 3 MHz period (0.33 µs), so the pulse passes almost undistorted.
\(h = c\Theta/2\pi\), and \(\Theta = 2\pi\) while the arrival circle is inside the aperture. For a disc the rectangle in \(h\) becomes \(p = \rho_0 c\, [v(t - z/c) - v(t - \sqrt{z^2 + a^2}/c)]\): two copies of the pulse, the second inverted. At z = 2 mm for a = 3 mm the edge wave lags by 1.04 µs, longer than the 3 MHz, 60 % pulse, so the two are separate; at 50 mm the lag is 0.058 µs and they merge into a derivative.
On axis the response runs from \(z/c\) to \(\sqrt{z^2 + a^2 + b^2}/c\), a duration \(\approx (a^2 + b^2)/(2zc)\). When that is much shorter than \(1/f_0\), \(h\) acts as a delta of area \(A/(2\pi R)\) and the pulse arrives unfiltered. For 4 × 6 mm at 3 MHz the crossover (\(\Delta t = 1/f_0\)) is at 12.7 mm; at 40 mm the pulse peak is 0.90 of the point-source value, at 15 mm only 0.44.
Each sub-element is replaced by a point-like response \(\Delta A/(2\pi R) \cdot \delta(t - R/c)\). That is exact only when the sub-element's own SIR is short compared with the pulse. For a 6 × 6 mm element at 3 MHz the on-axis error beyond 5 mm is 19 dB with 1.5 mm sub-elements, 3.1 dB with 0.75 mm, 0.37 dB with 0.375 mm and 0.08 dB with 0.19 mm (pulsed, 60 % bandwidth).
The disc formula \(\lvert p \rvert \propto \lvert \sin[(\pi/\lambda)(\sqrt{z^2 + a^2} - z)] \rvert\) has its last maximum where the path difference to the edge is \(\lambda/2\). For a square, the two edge pairs act as two Fresnel slits; their product peaks a little farther out. Computed from the closed-form on-axis SIR: 23.8 mm for a 6 mm square at 3 MHz (\(a^2/\lambda\) = 17.5 mm), \(1.36\, a^2/\lambda\); the Fresnel approximation, a product of two slit integrals, gives \(1.37\, a^2/\lambda\). With 60 % bandwidth the deep CW nulls fill in.
In a lossless plane wave the Fubini solution gives harmonic amplitudes \(B_n = 2J_n(n\sigma)/(n\sigma)\) with \(\sigma = x/\bar{x}\). In water (c = 1480 m/s, ρ = 1000 kg/m³, β = 3.5) at 3 MHz, \(\bar{x}\) is 98, 49 and 25 mm for 0.5, 1 and 2 MPa, so at 25 mm a 1 MPa wave already has a second harmonic 12.4 dB below the fundamental. Diagnostic foci reach these pressures, so measured harmonic content is real and a linear code will miss it.
Grid points grow as \((\text{domain}/\lambda)^3\) times \(\mathrm{PPW}^3\) and the number of time steps as \(\text{domain}/\lambda\) times \(\mathrm{PPW}/\mathrm{CFL}\). A 40 × 40 × 60 mm volume at 5 MHz is \(3.3 \times 10^6\, \lambda^3\): \(8.9 \times 10^{7}\) points at 3 PPW (about 3.5 GB for ten single-precision fields) but \(1.1 \times 10^{10}\) points and 440 GB at the 15 PPW a second-order scheme would need.
The FD2 dispersion relation \(\sin(\omega\Delta t/2) = \mathrm{CFL}\, \sin(k\Delta x/2)\) makes short wavelengths slow. At 10 PPW and CFL 0.3 the phase speed is 1.5 % low, 5.4° of phase per wavelength; at 20 PPW 0.37 %. FD4 at 8 PPW is within 0.06 %. In Demo 21.4 FD2 at 7.7 PPW widened a 2D focus by about 10 % and moved it 1–2 mm, an error the size of a real aberration.
FFT derivatives are exact for every resolvable wavenumber, so the spatial error is zero up to Nyquist. Plain leapfrog still distorts time (\(S = k\Delta x/2\) gives a fast phase speed: +0.23 % at 8 PPW and CFL 0.3). Multiplying the derivative by \(\operatorname{sinc}(c_{\mathrm{ref}} k \Delta t/2)\) makes \(\sin(\omega\Delta t/2) = \sin(c_{\mathrm{ref}} k \Delta t/2)\), i.e. exact propagation in a homogeneous medium for any step.
Matched broadband stacks often have a double-humped spectrum. In the figure two resonances at \(0.85 f_0\) and \(1.2 f_0\) (equal widths, the upper one 0.8 times as tall) give a peak at 4.33 MHz but a centre of 5.05 MHz for \(f_0\) = 5 MHz: the peak frequency misstates the band by 14 %. The arithmetic and geometric centres differ by \(f_c \cdot \mathrm{FBW}^2/8\) to first order (0.26 MHz at 5 MHz and 65 %), so state which one you report.
Echo spectra multiply the transmit and receive responses. With equal responses the two-way amplitude is \(\lvert H \rvert^2\), so its −6 dB points are the −3 dB points of \(\lvert H \rvert\). For a Gaussian band the one-way −6 dB bandwidth is \(\sqrt{2}\) times the two-way one: a probe quoted at 65 % pulse-echo radiates a one-way field with about 92 % at −6 dB. Mixing the two is a common datasheet error.
A Gaussian envelope of standard deviation \(\sigma_t\) has a Gaussian spectrum with \(\sigma_f = 1/(2\pi\sigma_t)\); each −6 dB width is \(2\sqrt{2 \ln 2}\) times its \(\sigma\), so the product is \(8 \ln 2/(2\pi) = 0.883\). At 5 MHz and 65 % the −6 dB envelope lasts 0.27 µs, an axial resolution of about 0.21 mm. A measured product well above 0.88 flags ringing or a non-Gaussian tail.
The measured echo is excitation × two-way insertion loss × reflection coefficient × water attenuation × diffraction loss. Steel in water reflects 0.937 in amplitude (−0.57 dB). At 5 MHz and a 30 mm stand-off the round-trip water loss is 0.33 dB, at 10 MHz 1.3 dB. Different reflectors, distances or receiver loads make raw loop sensitivities incomparable by a dB or more.
Near a jinc focus \(p \approx 1 - (\pi r/\lambda F\#)^2/8\); averaging \(r^2\) over a disc gives \(b^2/2\), hence the loss. For \(b = 0.2\, \lambda F\#\) it is 0.21 dB, for \(0.3\, \lambda F\#\) 0.48 dB, for \(0.5\, \lambda F\#\) 1.35 dB (exact). Intensity errors are double in dB. Harmonics and high-frequency probes have narrower beams and suffer more, which is why IEC 62127-1 ties the allowed size to wavelength and geometry.
A nonlinear wave steepens: its fronts are built from harmonics in phase, and a band-limited hydrophone both drops and phase-shifts them. For a plane Fubini waveform at σ = 0.8, a hydrophone with a first-order roll-off at \(3f_0\) reads \(p_+\) 9 % low and \(p_-\) 4 % low; at \(5f_0\) 3.2 % and 1.6 %; at \(10f_0\) under 1 %. Deconvolving with the complex sensitivity \(M_L(f)\) recovers most of it within the calibrated band.
Water scans have almost no loss, so the derating models tissue. At 3 MHz the factor is 0.9 dB per cm: −4.5 dB at 5 cm. For the focused piston in the figure (radius 6 mm, geometric focus 50 mm, 3 MPa peak in water at 35 mm because of focal shift) the derated maximum moves to 33 mm and MI falls from 1.73 (water value) to 1.22. The location of the derated maximum, not the water focus, defines MI.
The error weights are random with variance \(\sigma_a^2 w_n^2\), so their pattern has mean power \(\sigma_a^2 \sum w_n^2\) everywhere, against a peak \((\sum w_n)^2\). Uniform weights give \(\sigma_a^2/N\); Hann gives about \(1.5\, \sigma_a^2/N\). A dead element is an error of \(-w_n\). These floors are invisible under the −13 dB uniform sidelobes but dominate a low-sidelobe (apodized) design, which is why QA limits are tighter for apodized probes.
If each element drags its neighbours with amplitude \(\kappa\) in phase, the effective element is \([\kappa, 1, \kappa]\) wide and its pattern gains the factor \(1 + 2\kappa \cos(kp \sin\theta)\). At 0.3 mm pitch and 5 MHz, \(kp \sin\theta \approx \pi\) at 30°, so the factor falls from \(1 + 2\kappa\) to \(1 - 2\kappa\): −3.5 dB for −20 dB crosstalk, −6.4 dB for −15 dB. Phased arrays are therefore specified more tightly than linear arrays.
21.11Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Spatial impulse response | \(h = \int_S \delta(t - R/c)/(2\pi R)\, dS = c\Theta(t)/2\pi\) | Field II, FOCUS; any linear field |
| Pressure from SIR | \(p = \rho_0\, \partial v/\partial t * h\); \(P = j2\pi f \rho_0 V H\) | Pulsed and CW beam maps |
| Pulse-echo (point scatterer) | \(y \propto E_m * \partial^3 v/\partial t^3 * h_t * h_r\) | Image simulation |
| Disc on axis | \(p = \rho_0 c\, [v(t - z/c) - v(t - \sqrt{z^2 + a^2}/c)]\) | Analytic benchmark |
| Near-field length | disc \(a^2/\lambda\); square \(\approx 1.36\, a^2/\lambda\) | Where the beam starts to spread |
| Far-field condition (on axis) | \((a^2 + b^2)/(2zc) \ll 1/f_0\) | When point-source and Fraunhofer formulas hold |
| Leapfrog dispersion | \(\sin(\omega\Delta t/2) = \mathrm{CFL} \cdot S(k\Delta x)\) | Choosing PPW and CFL |
| FD2 phase-speed error | \(\approx -(k\Delta x)^2 (1 - \mathrm{CFL}^2)/24\) | 15–20 PPW for FD2 |
| k-space correction | \(jk \to jk\, \operatorname{sinc}(c_{\mathrm{ref}} k \Delta t/2)\) | k-Wave near 2–4 PPW |
| Shock distance | \(\bar{x} = \rho c^3/(\beta\omega p_0)\), \(B_n = 2J_n(n\sigma)/(n\sigma)\) | When linear models fail |
| Centre frequency, FBW | \((f_l + f_u)/2\); \((f_u - f_l)/f_c\) | Pulse-echo datasheets |
| Gaussian time-bandwidth | \(\tau_{-6} B_{-6} = 0.883\) | Sanity check of pulse length |
| Loop sensitivity | \(20 \log_{10}(V_{\mathrm{echo,pp}}/V_{\mathrm{tx,pp}})\) | Probe acceptance tests |
| Coupling from impedance | \(k_t^2 = (\pi/2)(f_r/f_a) \tan[(\pi/2)(f_a - f_r)/f_a]\) | Material and element QA |
| Hydrophone averaging | loss \(\approx (\pi b/\lambda F\#)^2/16\) | Field-scan corrections |
| Derating and MI | \(p_{\mathrm{r.3}} = p_r 10^{-0.3fz/20}\); \(\mathrm{MI} = p_{\mathrm{r.3}}/\sqrt{f}\) | IEC 62359, FDA tables |
| Element-error floor | \(\sigma_a^2 \sum w^2/(\sum w)^2\); \(D/(N - D)^2\) | Array QA limits |
| Crosstalk factor | \(1 + 2\kappa \cos(kp \sin\theta)\) | Steering loss in phased arrays |
21.12Interview and qualifying-exam questions
I1What is a spatial impulse response, and why is it the core of Field II? · screen
It is the velocity potential at a point when the aperture moves with an impulsive velocity: \(h(\mathbf{r}, t) = \int \delta(t - R/c)/(2\pi R)\, dS\). In a linear homogeneous medium every field quantity follows from it by convolution with the excitation (\(p = \rho_0 v' * h\)) and pulse-echo is \(h_t * h_r\). Field II splits the aperture into small rectangles and sums their far-field responses, so it never solves a wave equation. Follow-up: "What does it cost you?" No heterogeneity, no nonlinearity, no multiple scattering.
I2Sketch h(t) for a rectangular element at a point on axis and at a point beyond its edge. · onsite
On axis: zero until \(z/c\), a flat plateau at \(c\) until the circle reaches the nearer edge (\(\sqrt{z^2 + a^2}/c\)), a decay through the farther edge, and zero after the corners at \(\sqrt{z^2 + a^2 + b^2}/c\). Beyond the edge: starts when the circle touches the nearest edge, never reaches \(c\) (only part of the circle is inside), shows knees at the corner arrivals, and lasts much longer. For 4 × 6 mm at z = 15 mm: on axis 0.28 µs at height \(c\); at x = 6 mm, 1.13 µs with peak 0.20c. Follow-up: "What does a long h do to the pulse?" It low-pass filters it, so off-axis pulses are longer and weaker.
I3When would you use k-Wave instead of Field II for a sound-speed estimation study? · research
Whenever the estimator's behaviour depends on the heterogeneity: layered fat and muscle, refraction at oblique interfaces, aberration-induced coherence loss, reverberation clutter. Field II with a global c only tests the homogeneous mismatch problem. I would validate the k-Wave setup in homogeneous water against Field II first, choose PPW from a convergence test, and smooth sharp interfaces. Follow-up: "Cost?" A 2D study at 5 MHz is cheap; 3D needs a GPU and tens of GB at 3 PPW.
I4Your FDTD simulation of a homogeneous medium shows the focus 1.5 mm shallower than expected. What do you check? · onsite
Numerical dispersion first: at 7–8 PPW a second-order scheme has phase-speed errors of 2–5 % across the band and they depend on direction, which distorts focusing (Demo 21.4). Then CFL, source discretization, absorbing boundary reflections, and whether the delays use the same c as the grid. Refine \(\Delta x\) by 2 and see if the focus moves; switch to fourth order or pseudospectral. Follow-up: "What is the convergence rate?" \(\Delta x^2\) for FD2, so halving the grid should cut the error by about 4.
I5Define centre frequency and fractional bandwidth from a pulse-echo test. What mistakes are common? · screen
Take the amplitude spectrum of the gated echo, find the −6 dB points \(f_l\), \(f_u\); \(f_c = (f_l + f_u)/2\), \(\mathrm{FBW} = (f_u - f_l)/f_c\). Mistakes: reporting the peak as the centre (wrong by 10–20 % for double-humped spectra), mixing −3 dB one-way with −6 dB two-way, gates that include reverberations (ripple) or cut the echo (broadening), power versus amplitude dB, and not stating arithmetic versus geometric centre. Follow-up: "How long is a 5 MHz, 65 % pulse?" About 0.27 µs at −6 dB for a Gaussian.
I6What is loop sensitivity, and how would you compare it between two probes measured in different labs? · onsite
\(20 \log_{10}(V_{\mathrm{echo,pp}}/V_{\mathrm{tx,pp}})\) for a flat reflector. To compare, normalize the reflector (steel −0.57 dB, quartz differs), the water path (\(2.2 \times 10^{-4}\) dB/(mm·MHz²) round trip), diffraction to the reflector distance, the pulser source impedance and the receiver load (50 Ω versus high impedance changes the voltage), and the cable. Follow-up: "Is loop sensitivity a transducer property?" No, it is a property of transducer + cable + electronics + target geometry.
I7How large may a hydrophone be for a focused 10 MHz beam at F# 2? · research
\(\lambda F\# = 0.31\) mm. For under 0.5 dB peak loss the radius should be below about \(0.3\, \lambda F\#\): \(b \lesssim 0.09\) mm, a diameter under about 0.2 mm, and harmonics need smaller still. The loss law is \((\pi b/\lambda F\#)^2/16\). Then check the hydrophone's directivity, its effective (not geometric) size at the frequency, and its bandwidth for the harmonic content. IEC 62127-1 gives the formal size criterion and corrections. Follow-up: "And for \(I_{\mathrm{SPTA}}\)?" The intensity error is twice the pressure error in dB.
I8Walk me through how MI is computed for a 510(k). · onsite
Hydrophone scans in water under IEC 62127-1 conditions; deconvolve or correct for sensitivity, bandwidth and spatial averaging; derate the peak rarefactional pressure by 0.3 dB/cm/MHz; find the depth defined by IEC 62359 (via the derated pulse-intensity integral); \(\mathrm{MI} = p_{\mathrm{r.3}}/\sqrt{f_{\mathrm{awf}}}\) with p in MPa and f in MHz; search operating modes for the global maximum; report against the FDA Track 3 limit MI ≤ 1.9 with uncertainty. Follow-up: "Why \(\sqrt{f}\)?" The cavitation threshold data scale roughly as \(p \propto \sqrt{f}\), so MI approximates the cavitation likelihood.
I9A probe returns from the field with two dead elements. Is it still usable? · screen
Quantify: for 64 active elements, uniform weights, the main lobe drops \(20 \log(62/64)\) = 0.28 dB and a far-sidelobe floor appears near \(10 \log(2/62^2)\) = −32.8 dB; with an apodized aperture the floor is relative to a much lower sidelobe level and is more visible. Two adjacent dead elements produce structured lobes rather than a floor. Whether that fails depends on the application's contrast requirement and the QA limit; many services set a maximum number of dead elements and forbid adjacent ones. Follow-up: "How do you detect them?" Per-element pulse-echo or impedance; a phantom image shows vertical dropout streaks in linear arrays.
I10How do you turn an impedance measurement into a KLM model? · research
Start with the bare plate in air: \(f_r\) and \(f_a\) give \(k_t\) through the IEEE 176 relation, the capacitance well below resonance times \(1 - k_t^2\) (or the capacitance well above resonance) gives \(\varepsilon^S\), the \(Q\) of the peaks gives the losses; \(f_a\) fixes the plate's sound speed for its thickness. Then add layers and backing and fit their impedances and thicknesses to the finished element's impedance in air and water, and finally fit the electrical network to the pulse-echo spectrum. Peaks the 1D model cannot produce are lateral modes: hand those to FE. Follow-up: "Example?" \(f_r\) = 3.00, \(f_a\) = 3.45 MHz gives \(k_t\) = 0.53.
21.13MATLAB lab
Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run.
m21_spatial_impulse_response.mExact SIR of a rectangle on and off axis, Field II-style subdivision, pressure pulses and on-axis field · change w, H, the field points and the sub-element count · runtime ≈ 4 s
%% m21_spatial_impulse_response.m
% Spatial impulse response (SIR) of a baffled rectangular element, exact and Field II-style.
% Everything here is computed by this script; Field II itself is not called (the subdivision
% only reimplements its far-field sub-element idea).
% h(r,t) = (c/2pi) * Theta(t), Theta = angle of the circle of radius sqrt((ct)^2 - z^2)
% (centred on the projection of r) that lies inside the aperture. p = rho * dv/dt * h.
% What it shows
% 1) exact SIR on axis, inside the edge and beyond the edge; check against the closed form
% 2) integral of h equals the surface integral of 1/(2 pi R)
% 3) Field II-style subdivision (point-like sub-elements) converging to the exact SIR
% 4) on-axis CW field of a 6 mm square: subdivided sum versus the exact SIR
% What to change: element size (w, H), field points, f0, sub-element counts.
% Expected output: plateau h/c = 1.000 on axis, off-axis duration 1.129 us and peak 0.204c,
% last on-axis CW maximum of the 6 mm square near 24 mm, errors ~22, 4.4, 0.6, 0.13 dB.
% Units: mm, us, MHz; c in mm/us. Base MATLAB / Octave.
clear; close all; rng(1);
c = 1.54; % mm/us (1540 m/s)
w = 4; H = 6; % element width (x) and height (y), mm
a = w/2; b = H/2;
fs = 200; % MHz, SIR sampling
nAng = 1440; % angles on the arrival circle
phi = ((0:nAng-1) + 0.5) * 2*pi/nAng;
% exact SIR by testing points on the arrival circle (vectorized over angles)
sirRect = @(x0, y0, z0, t) c * mean( ...
(abs(x0 + sqrt(max((c*t(:)).^2 - z0^2, 0)) * cos(phi)) < a) & ...
(abs(y0 + sqrt(max((c*t(:)).^2 - z0^2, 0)) * sin(phi)) < b) & ...
repmat((c*t(:)).^2 > z0^2, 1, nAng), 2).'; % returns a row vector
%% 1) on axis, inside the edge, beyond the edge (z = 15 mm)
z0 = 15;
pts = [0 0; 1.5 0; 6 0];
names = {'on axis', 'x = 1.5 mm (inside)', 'x = 6 mm (beyond edge)'};
f0 = 3; B = 0.6; % pulse: 3 MHz, 60 % -6 dB bandwidth
sf = B*f0/(2*sqrt(2*log(2))); st = 1/(2*pi*sf);
dv = @(t) exp(-t.^2/(2*st^2)) .* (2*pi*f0*cos(2*pi*f0*t) - t/st^2 .* sin(2*pi*f0*t));
tk = (-ceil(4*st*fs):ceil(4*st*fs))/fs; % pulse kernel time axis
figure(1); clf;
for k = 1:3
x0 = pts(k,1); y0 = pts(k,2);
Rmin = sqrt(max(abs(x0)-a,0)^2 + max(abs(y0)-b,0)^2 + z0^2);
Rmax = sqrt((abs(x0)+a)^2 + (abs(y0)+b)^2 + z0^2);
t = (floor(Rmin/c*fs)-2 : ceil(Rmax/c*fs)+2)/fs;
h = sirRect(x0, y0, z0, t);
corners = sort(sqrt((x0-[-a a -a a]).^2 + (y0-[-b -b b b]).^2 + z0^2)/c);
fprintf('%-24s t0 = %.3f us, t1 = %.3f us, duration = %.3f us, peak h/c = %.3f\n', ...
names{k}, Rmin/c, Rmax/c, (Rmax-Rmin)/c, max(h)/c);
fprintf('%-24s corner arrivals: %s us\n', '', sprintf('%.3f ', corners));
p = conv(h, dv(tk)) / fs; % rho = 1
tp = t(1) + tk(1) + (0:numel(p)-1)/fs;
pff = w*H/(2*pi*sqrt(x0^2+y0^2+z0^2)) * max(abs(dv(tk)));
fprintf('%-24s pressure peak / point-source estimate = %.3f\n', '', max(abs(p))/pff);
subplot(2,1,1); plot(t - z0/c, h/c, 'LineWidth', 1.5); hold on;
subplot(2,1,2); plot(tp - z0/c, p/pff, 'LineWidth', 1.2); hold on;
if k == 1
% closed form on axis: Theta = 2pi - 4acos(a/s)[s>a] - 4acos(b/s)[s>b]
s = sqrt(max((c*t).^2 - z0^2, 0));
th = 2*pi*ones(size(s));
th(s > a) = th(s > a) - 4*acos(a./s(s > a));
th(s > b) = th(s > b) - 4*acos(b./s(s > b));
th(s >= sqrt(a^2+b^2) | (c*t).^2 <= z0^2) = 0;
fprintf('%-24s max |numeric - closed form| / c = %.1e\n', '', max(abs(h - c*th/(2*pi)))/c);
end
end
subplot(2,1,1); ylim([0 1.15]); xlabel('t - z/c (us)'); ylabel('h / c'); title('Spatial impulse response, 4 x 6 mm, z = 15 mm'); legend(names); grid on;
subplot(2,1,2); xlabel('t - z/c (us)'); ylabel('p / point-source peak'); title('Pressure p = \rho v'' * h, 3 MHz 60 % pulse'); grid on;
%% 2) integral check: int h dt = int dS / (2 pi R)
x0 = 5; y0 = 1; z0 = 8;
Rmin = sqrt(max(abs(x0)-a,0)^2 + max(abs(y0)-b,0)^2 + z0^2); Rmax = sqrt((abs(x0)+a)^2 + (abs(y0)+b)^2 + z0^2);
t = (floor(Rmin/c*800):ceil(Rmax/c*800))/800;
I1 = sum(sirRect(x0, y0, z0, t)) / 800;
ng = 400; [XS, YS] = meshgrid(-a + ((1:ng)-0.5)*w/ng, -b + ((1:ng)-0.5)*H/ng);
I2 = sum(sum((w/ng)*(H/ng) ./ (2*pi*sqrt((x0-XS).^2 + (y0-YS).^2 + z0^2))));
fprintf('\nIntegral check at (5,1,8) mm: int h dt = %.5f mm, surface integral = %.5f mm\n', I1, I2);
%% 3) Field II-style subdivision of the SIR at the off-axis point (6, 0, 15)
x0 = 6; y0 = 0; z0 = 15; fsb = 200;
Rmin = sqrt((x0-a)^2 + z0^2); Rmax = sqrt((x0+a)^2 + b^2 + z0^2);
t = (floor(Rmin/c*fsb)-4 : ceil(Rmax/c*fsb)+4)/fsb;
hx = sirRect(x0, y0, z0, t);
pex = conv(hx, dv(tk))/fsb;
fprintf('\nSubdivision at (6, 0, 15) mm, pressure-pulse error vs exact:\n');
for n = [2 4 8 16 32]
[XS, YS] = meshgrid(-a + ((1:n)-0.5)*w/n, -b + ((1:n)-0.5)*H/n);
R = sqrt((x0-XS(:)).^2 + (y0-YS(:)).^2 + z0^2); amp = (w/n)*(H/n) ./ (2*pi*R);
u = (R/c - t(1))*fsb; i0 = floor(u); fr = u - i0; % linear split of each delta
hs = accumarray([i0+1; i0+2], [amp.*(1-fr); amp.*fr], [numel(t)+2, 1])' * fsb;
hs = hs(1:numel(t));
ps = conv(hs, dv(tk))/fsb;
fprintf(' %2d x %2d sub-elements (%.2f x %.2f mm): rms error %.1f %% of peak\n', n, n, w/n, H/n, 100*sqrt(mean((ps-pex).^2))/max(abs(pex)));
end
%% 4) on-axis CW field of a 6 mm square at 3 MHz: exact SIR vs subdivision
ws = 6; as = 3; f = 3; k = 2*pi*f/c; lam = c/f;
zs = 2:0.25:60;
Pex = zeros(size(zs));
for iz = 1:numel(zs)
z = zs(iz); t0 = z/c; t1 = sqrt(z^2 + 2*as^2)/c; nt = max(400, ceil((t1-t0)*2000));
tt = t0 + ((1:nt)-0.5)*(t1-t0)/nt; s = sqrt((c*tt).^2 - z^2);
th = 2*pi - 8*acos(min(as./s, 1)); % square: a = b
th(s >= sqrt(2)*as) = 0;
Pex(iz) = abs(sum(c*th/(2*pi) .* exp(-1i*2*pi*f*tt)) * (t1-t0)/nt);
end
mx = false(size(Pex)); WN = 24;
for i = WN+1:numel(Pex)-WN
mx(i) = Pex(i) == max(Pex(i-WN:i+WN));
end
iz = find(mx, 1, 'last');
fprintf('\n6 mm square, 3 MHz: last on-axis CW maximum at z = %.2f mm (a^2/lambda = %.2f mm, ratio %.2f)\n', zs(iz), as^2/lam, zs(iz)/(as^2/lam));
figure(2); clf; plot(zs, 20*log10(Pex/max(Pex)), 'k', 'LineWidth', 1.8); hold on;
for n = [1 8 16 32]
[XS, YS] = meshgrid(-as + ((1:n)-0.5)*ws/n, -as + ((1:n)-0.5)*ws/n); dA = (ws/n)^2;
R = sqrt(XS(:).^2 + YS(:).^2 + zs.^2); % (n^2) x numel(zs)
Pap = abs(sum(dA ./ (2*pi*R) .* exp(-1i*k*R), 1));
e = max(abs(20*log10(Pap(zs >= 5) ./ Pex(zs >= 5))));
fprintf(' %2d x %2d sub-elements (%.3f mm = %.2f lambda): max on-axis error (z >= 5 mm) = %.2f dB\n', n, n, ws/n, ws/n/lam, e);
plot(zs, 20*log10(Pap/max(Pex)), '--', 'LineWidth', 1.1);
end
xlabel('z (mm)'); ylabel('|p| on axis (dB re max)'); ylim([-30 3]); grid on;
title('6 mm square, 3 MHz CW: exact SIR (black) vs subdivided sums');
legend('exact', '1x1', '8x8', '16x16', '32x32', 'Location', 'southeast');
What to look for: the on-axis plateau prints \(h/c\) = 1.000 and matches the closed form to a few parts in 10³ (the error sits at the knee, a sampling effect); the off-axis point (x = 6 mm) lasts 1.13 µs with peak 0.20c; \(\int h\, dt\) agrees with the direct surface integral; the subdivided on-axis CW field approaches the exact one as sub-elements shrink (errors of order 20, 4, 0.6 and 0.1 dB for 1, 8, 16 and 32 per side); the last on-axis CW maximum of the 6 mm square is at 23.75 mm (\(1.35\, a^2/\lambda\) on the 0.25 mm grid); the subdivided off-axis pressure pulse converges as 8 % → 1.6 % → 0.4 % rms error for 8, 16, 32 sub-elements per side.
m21_pulse_echo_analysis.mSynthetic measured echo, spectrum and envelope analysis, loop sensitivity, gate-length sweep · change fc, bandwidth, reverberation and gate · runtime ≈ 1 s
%% m21_pulse_echo_analysis.m
% Pulse-echo test analysis on a synthetic "measured" echo from a flat steel reflector.
% What it shows
% centre frequency (arithmetic and geometric), -6 dB bandwidth, -6 / -20 dB pulse length,
% loop sensitivity and its correction for the reflector and water path, and how the gate
% length changes the answer when a lens reverberation is present.
% What to change: fc, bw (two-way -6 dB FBW), reverberation level/delay, SNR, gate length.
% Expected output (defaults): fc ~ 5.0 MHz, FBW ~ 65 %, loop sensitivity ~ -20.9 dB for a true
% two-way insertion loss of -20 dB (steel -0.57 dB, water -0.33 dB); FBW jumps once the gate
% includes the reverberation (-14 dB at 1.6 us) and the spectrum ripples with period 1/1.6 us.
% Units: us, MHz, V. Base MATLAB / Octave.
clear; close all; rng(1);
fs = 250; N = 2000; t = (0:N-1)/fs; % 250 MS/s scope, 8 us record
fc = 5; bw = 0.65; % transducer: centre frequency (MHz), two-way -6 dB FBW
rv_dB = -14; dl = 1.6; % lens reverberation level (dB) and delay (us)
snr_dB = 40; IL = -20; % echo peak / rms noise; true two-way insertion loss (dB)
Vtx = 100; z = 30; % excitation (Vpp), reflector distance (mm)
Rst = (45.4 - 1.48)/(45.4 + 1.48); % steel / water amplitude reflection
aw = 2.2e-4; % water attenuation, dB/(mm MHz^2)
TE = 1.0; % echo arrival within the record (us)
% causal gamma-tone echo, n = 3: -6 dB FBW = 0.2049/(fc*tau)
n = 3; tau = 0.2049/(bw*fc);
g = @(tt) (tt > 0) .* (max(tt,0)/tau).^n .* exp(-max(tt,0)/tau) .* cos(2*pi*fc*tt);
e0 = g(t - TE);
Vpp = Vtx * 10^(IL/20) * Rst * 10^(-2*z*aw*fc^2/20);
x = (e0 - 10^(rv_dB/20) * g(t - TE - dl)) * Vpp / (max(e0) - min(e0));
x = x + (Vpp/2) * 10^(-snr_dB/20) * randn(size(x));
% analysis function (anonymous helpers only)
nfft = 16384;
tukey = @(L, nt) [0.5-0.5*cos(pi*(0:nt-1)/nt), ones(1, L-2*nt), 0.5-0.5*cos(pi*(nt-1:-1:0)/nt)];
gates = [1.2 0.6 0.8 1.0 1.6 2.0 2.5 3.0];
fprintf('gate (us) fpk fl fu fc(ar) fc(geo) FBW(%%) d-6(us) d-20(us)\n');
for ig = 1:numel(gates)
gl = gates(ig);
i0 = round((TE - 0.25)*fs) + 1; L = round(gl*fs); nt = max(2, round(min(0.08*L, 0.1*fs)));
seg = x(i0:i0+L-1) .* tukey(L, nt);
X = fft(seg, nfft); mag = abs(X(1:nfft/2+1)); df = fs/nfft; f = (0:nfft/2)*df;
[mx, ip] = max(mag(2:end)); ip = ip + 1; thr = mx/2;
kl = ip; while kl > 2 && mag(kl) >= thr, kl = kl - 1; end
ku = ip; while ku < numel(mag) && mag(ku) >= thr, ku = ku + 1; end
fl = f(kl) + (thr - mag(kl))/(mag(kl+1) - mag(kl))*df;
fu = f(ku-1) + (mag(ku-1) - thr)/(mag(ku-1) - mag(ku))*df;
fcA = (fl + fu)/2; fcG = sqrt(fl*fu); fbw = (fu - fl)/fcA;
% envelope by the FFT analytic signal
Xa = X; Xa(2:nfft/2) = 2*Xa(2:nfft/2); Xa(nfft/2+2:end) = 0; xa = ifft(Xa); env = abs(xa(1:L));
[em, ie] = max(env);
lev = @(dbv) em*10^(dbv/20);
ia = ie; while ia > 1 && env(ia) >= lev(-6), ia = ia - 1; end
ib = ie; while ib < L && env(ib) >= lev(-6), ib = ib + 1; end
d6 = (ib - ia)/fs;
ia = ie; while ia > 1 && env(ia) >= lev(-20), ia = ia - 1; end
ib = ie; while ib < L && env(ib) >= lev(-20), ib = ib + 1; end
d20 = (ib - ia)/fs;
fprintf('%6.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.1f %5.3f %5.3f\n', gl, f(ip), fl, fu, fcA, fcG, 100*fbw, d6, d20);
if ig == 1
LS = 20*log10((max(seg) - min(seg))/Vtx);
ILc = LS - 20*log10(Rst) + 2*z*aw*fcA^2;
d6P = d6; segP = seg; envP = env; tP = t(i0:i0+L-1); fP = f; magP = mag; fl1 = fl; fu1 = fu; fc1 = fcA; fbw1 = fbw;
end
end
fprintf('\nGaussian time-bandwidth check: 0.883/(FBW*fc) = %.3f us; the causal gamma-tone echo measures %.3f us\n', 0.8825/(fbw1*fc1), d6P);
fprintf('Loop sensitivity (gate 1.2 us) = %.2f dB re %d Vpp\n', LS, Vtx);
fprintf('Corrected for steel (%.2f dB) and water (%.2f dB): %.2f dB (true %.1f dB)\n', 20*log10(Rst), -2*z*aw*fc1^2, ILc, IL);
fprintf('Arithmetic - geometric centre = %.3f MHz; fc*FBW^2/8 = %.3f MHz\n', fc1 - sqrt(fl1*fu1), fc1*fbw1^2/8);
figure(1); clf;
subplot(2,1,1); plot(t, x, 'b'); hold on; plot(tP, envP, 'Color', [0 0.5 0.4], 'LineWidth', 1.5);
yl = [-0.7 0.7]*Vpp; plot([tP(1) tP(1)], yl, '--', 'Color', [0.8 0.55 0]); plot([tP(end) tP(end)], yl, '--', 'Color', [0.8 0.55 0]); ylim(yl);
xlabel('t (us)'); ylabel('V'); title('Measured-like echo with gate (1.2 us) and envelope'); xlim([0 5]); grid on;
subplot(2,1,2); plot(fP, 20*log10(magP/max(magP)), 'b', 'LineWidth', 1.5); hold on;
plot([0 2*fc], [-6.02 -6.02], 'r--'); plot([fl1 fu1], [-6.02 -6.02], 'ro', 'MarkerFaceColor', 'r');
xlim([0 2*fc]); ylim([-40 2]); xlabel('f (MHz)'); ylabel('dB re peak'); grid on;
title(sprintf('Spectrum: fc = %.2f MHz, FBW = %.1f %%', fc1, 100*fbw1));
What to look for: with the default 5 MHz, 65 % echo the analysis returns \(f_c\) = 5.00 MHz and FBW = 64.8 % for the 1.2 µs gate; a 0.6 µs gate cuts the echo and reports 95 %; gates of 2.5 µs and longer include the −14 dB reverberation at 1.6 µs and the −6 dB band collapses to 47 %; the loop sensitivity prints about −20.9 dB for a true insertion loss of −20 dB (reflector −0.57 dB, water −0.33 dB).
m21_hydrophone_averaging.mSpatial averaging over a disc at a jinc focus, the small-radius law, harmonic beams, and band-limited measurement of a Fubini waveform · change frequency, F#, diameter, roll-off · runtime ≈ 2 s
%% m21_hydrophone_averaging.m
% Hydrophone spatial averaging at a focus and band-limited measurement of a steepened wave.
% What it shows
% 1) a hydrophone of radius b reads the disc average of p(r) = jinc(r/(lambda F#));
% loss ~ (pi b / lambda F#)^2 / 16 for small b; harmonics (~p^2) suffer more
% 2) a worked case (5 MHz, F# 2, 0.4 mm) and the same hydrophone at 15 MHz
% 3) a plane Fubini waveform (sigma = 0.8) read through a first-order hydrophone roll-off
% What to change: frequency, F#, diameter, sigma, roll-off frequencies.
% Expected output: -0.215 dB at b/lambdaF# = 0.2, -1.35 dB at 0.5; 0.4 mm at 5 MHz F#2: -0.57 dB;
% at 15 MHz about -5.2 dB; 3f0 roll-off: p+ about 9 % low, p- about 4 % low.
% Units: mm, MHz. Base MATLAB / Octave (besselj is part of base MATLAB).
clear; close all; rng(1);
c = 1.54;
jinc = @(x) (abs(x) < 1e-9) + (abs(x) >= 1e-9) .* 2 .* besselj(1, pi*x) ./ (pi*x + (abs(x) < 1e-9));
nr = 40; nt = 64; % polar quadrature over the disc
rq = sqrt(((1:nr) - 0.5)/nr); tq = ((1:nt) - 0.5)*2*pi/nt;
[RQ, TQ] = meshgrid(rq, tq);
diskAvg = @(prof, x0, bb) mean(mean(prof(sqrt((x0 + bb*RQ.*cos(TQ)).^2 + (bb*RQ.*sin(TQ)).^2))));
%% 1) loss table versus normalized radius
fprintf('b/(lambda F#) exact (dB) small-b law (dB) 2nd harmonic (dB)\n');
bs = [0.05 0.1 0.2 0.3 0.4 0.5];
for bb = bs
e1 = 20*log10(diskAvg(jinc, 0, bb));
ap = 20*log10(1 - (pi*bb)^2/16);
e2 = 20*log10(diskAvg(@(r) jinc(r).^2, 0, bb));
fprintf(' %4.2f %7.3f %7.3f %7.3f\n', bb, e1, ap, e2);
end
%% 2) worked cases
cases = [5 2 0.4; 15 2 0.4; 7.5 1.5 0.2; 10 2 0.2];
fprintf('\nf (MHz) F# d (mm) lambdaF# (mm) b/lambdaF# loss (dB) pressure intensity\n');
for k = 1:size(cases, 1)
f = cases(k,1); Fn = cases(k,2); d = cases(k,3); lF = c/f*Fn; bb = d/2/lF;
m = diskAvg(jinc, 0, bb);
fprintf('%5.1f %4.1f %5.2f %6.3f %5.3f %6.2f %5.1f %% %5.1f %%\n', f, Fn, d, lF, bb, 20*log10(m), 100*(1-m), 100*(1-m^2));
end
% measured profile for the first case
f = 5; Fn = 2; d = 0.4; lF = c/f*Fn; bb = d/2/lF;
u = linspace(-2, 2, 201); tr = jinc(abs(u)); ms = zeros(size(u));
for i = 1:numel(u), ms(i) = diskAvg(jinc, u(i), bb); end
ml = 101:(100 + find(tr(101:end) < 0, 1)); % main lobe (positive half)
hw = @(p) 2*interp1(p(ml) - p(101)/2, u(ml), 0); % -6 dB full width
fprintf('\n-6 dB width true %.3f mm, measured %.3f mm (1.41 lambda F# = %.3f mm)\n', hw(tr)*lF, hw(ms)*lF, 1.4102*lF);
figure(1); clf;
subplot(2,1,1); plot(u*lF, tr, 'k', u*lF, ms, 'b', 'LineWidth', 1.5); grid on;
xlabel('x (mm)'); ylabel('p (norm.)'); title('Focal profile, 5 MHz, F# 2: true (black) and 0.4 mm hydrophone (blue)');
bb2 = linspace(0, 0.6, 31); L1 = zeros(size(bb2)); L2 = L1;
for i = 1:numel(bb2), L1(i) = 20*log10(diskAvg(jinc, 0, bb2(i))); L2(i) = 20*log10(diskAvg(@(r) jinc(r).^2, 0, bb2(i))); end
subplot(2,1,2); plot(bb2, L1, 'b', bb2, 20*log10(1 - (pi*bb2).^2/16), 'k--', bb2, L2, 'r', 'LineWidth', 1.5); grid on;
xlabel('b / (\lambda F#)'); ylabel('peak loss (dB)'); legend('fundamental', 'small-b law', '2nd harmonic', 'Location', 'southwest');
%% 3) band-limited measurement of a Fubini waveform
sig = 0.8; NH = 60; tt = (0:4095)/4096; % one period
Bn = zeros(1, NH); for n = 1:NH, Bn(n) = 2*besselj(n, n*sig)/(n*sig); end
pt = Bn * sin(2*pi*(1:NH)'*tt);
fprintf('\nFubini sigma = %.2f: B2/B1 = %.3f (%.1f dB)\n', sig, Bn(2)/Bn(1), 20*log10(Bn(2)/Bn(1)));
figure(2); clf; plot(tt, pt, 'k', 'LineWidth', 1.6); hold on;
for fcut = [3 5 10]
nn = (1:NH)'; mag = 1./sqrt(1 + (nn/fcut).^2); ph = -atan(nn/fcut);
pm = (Bn' .* mag)' * sin(2*pi*nn*tt + repmat(ph, 1, numel(tt)));
fprintf('roll-off %2d f0: p+ %5.1f %% low, p- %5.1f %% low\n', fcut, 100*(1 - max(pm)/max(pt)), 100*(1 - min(pm)/min(pt)));
plot(tt, pm, 'LineWidth', 1.1);
end
xlabel('t / period'); ylabel('p / p_0'); grid on; title('Fubini waveform, \sigma = 0.8: true and band-limited hydrophone');
legend('true', '3 f_0', '5 f_0', '10 f_0', 'Location', 'southwest');
What to look for: the table of peak loss versus \(b/\lambda F\#\) matches \((\pi b/\lambda F\#)^2/16\) for small radii (−0.215 dB at 0.2, −1.35 dB exact at 0.5); the 5 MHz, F# 2, 0.4 mm case loses about 0.6 dB; the Fubini waveform at σ = 0.8 read through a first-order \(3f_0\) roll-off under-reads \(p_+\) by about 9 % and \(p_-\) by about 4 %.
m21_numerical_dispersion.m1D time-stepping with FD2, FD4 and k-space pseudospectral schemes over 100 wavelengths, measured delay versus the dispersion relation · change PPW and CFL · runtime ≈ 2 s
%% m21_numerical_dispersion.m
% 1D time stepping with three schemes over 100 wavelengths, compared with their dispersion relations.
% staggered leapfrog: sin(w dt/2) = CFL * S(k dx)
% FD2 S = sin(q/2); FD4 S = 9/8 sin(q/2) - 1/24 sin(3q/2); k-space PSTD: exact (w = c k)
% What it shows: the measured phase lag at f0 after 100 lambda equals the dispersion-relation
% prediction; FD2 at 8 PPW is more than two periods late, FD4 about 0.06 cycle (~20 deg) early, k-space exact.
% What to change: PPW, CFL, distance.
% Expected output (8 PPW, CFL 0.3): FD2 -2.34 cycles, FD4 +0.06, k-space 0.00 (predicted = measured).
% Units: wavelength lambda0 = 1, c = 1, rho = 1. Base MATLAB / Octave.
% These are small 1D solvers written for this course, not k-Wave; the k-space scheme uses the
% same correction idea (Tabei, Mast & Waag 2002) on a periodic 1D grid.
clear; close all; rng(1);
PPW = 8; CFL = 0.3; Ldist = 100; % points per wavelength at f0, Courant number, distance (lambda)
Ndom = 128; % domain length in wavelengths (periodic)
N = Ndom*PPW; dx = 1/PPW; dt = CFL*dx; nsteps = round(Ldist/dt); T = nsteps*dt;
x = (0:N-1)*dx; % pressure nodes
xv = x + dx/2; % velocity nodes (staggered)
f0 = 1; B = 0.6; sf = B*f0/(2*sqrt(2*log(2))); st = 1/(2*pi*sf);
x0 = 20; % initial pulse centre
pulse = @(xx) exp(-(xx - x0).^2/(2*st^2)) .* sin(2*pi*f0*(xx - x0));
k = 2*pi*[0:N/2, -N/2+1:-1]/(N*dx); % FFT wavenumbers
sinc_ = @(z) (abs(z) < 1e-12) + (abs(z) >= 1e-12) .* sin(z) ./ (z + (abs(z) < 1e-12));
q0 = 2*pi/PPW; k0bin = Ndom + 1; % f0 sits exactly on bin Ndom (domain = 128 lambda)
schemes = {'FD2', 'FD4', 'k-space'};
Sfun = {@(q) sin(q/2), @(q) 9/8*sin(q/2) - 1/24*sin(3*q/2), []};
figure(1); clf;
for is = 1:3
p = pulse(x);
v = pulse(xv + dt/2); % right-going p(x - ct): v = p/(rho c) at t = -dt/2
tic;
switch is
case 1
for n = 1:nsteps
v = v - CFL*(circshift(p, [0 -1]) - p);
p = p - CFL*(v - circshift(v, [0 1]));
end
case 2
A1 = 9/8; A2 = -1/24;
for n = 1:nsteps
v = v - CFL*(A1*(circshift(p, [0 -1]) - p) + A2*(circshift(p, [0 -2]) - circshift(p, [0 1])));
p = p - CFL*(A1*(v - circshift(v, [0 1])) + A2*(circshift(v, [0 -1]) - circshift(v, [0 2])));
end
case 3
kap = sinc_(k*dt/2); % k-space correction, c_ref = 1
sh = exp(1i*k*dx/2); % staggering shift
for n = 1:nsteps
v = v - dt*real(ifft(1i*k.*kap.*sh .* fft(p)));
p = p - dt*real(ifft(1i*k.*kap.*conj(sh) .* fft(v)));
end
end
el = toc;
% predicted phase error at f0 (cycles): (w_num - c k0) T / 2pi
if is < 3
wn = 2*asin(CFL*Sfun{is}(q0))/dt;
else
wn = 2*pi*f0;
end
pred = (wn - 2*pi*f0)*T/(2*pi);
% measured: phase of bin f0 relative to the exact translation
P0 = fft(pulse(x)); P1 = fft(p);
expct = P0(k0bin) * exp(-1i*2*pi*f0*T);
meas = angle(P1(k0bin)/expct)/(2*pi); % cycles, wrapped to (-0.5, 0.5]
wrapP = pred - round(pred);
fprintf('%-8s predicted %+7.3f cycles (wrapped %+6.3f), measured (wrapped) %+6.3f, travel-time error %+6.3f %%, %.1f s\n', ...
schemes{is}, pred, wrapP, -meas, -100*pred/Ldist, el);
xs = mod(x - x0 - T + Ndom/2, Ndom) - Ndom/2; % moving frame, centred
[xs, io] = sort(xs);
subplot(3,1,is); plot(xs, p(io), 'LineWidth', 1.2); hold on; plot(xs, pulse(xs + x0), 'k--');
xlim([-8 4]); grid on; ylabel('p'); title(sprintf('%s after %d \\lambda at %g PPW, CFL %.2f', schemes{is}, Ldist, PPW, CFL));
end
xlabel('(x - ct)/\lambda_0');
% dispersion curves
ppw = linspace(3, 20, 200); qq = 2*pi./ppw;
figure(2); clf;
plot(ppw, 100*(2*asin(CFL*sin(qq/2))./(CFL*qq) - 1), 'r', ppw, 100*(2*asin(CFL*(9/8*sin(qq/2) - 1/24*sin(3*qq/2)))./(CFL*qq) - 1), 'Color', [0.7 0.45 0], 'LineWidth', 1.6); hold on;
plot(ppw, 100*(2*asin(min(1, CFL*qq/2))./(CFL*qq) - 1), 'b', ppw, 0*ppw, 'g--', 'LineWidth', 1.6);
ylim([-8 3]); grid on; xlabel('points per wavelength'); ylabel('phase-speed error (%)');
legend('FD2', 'FD4', 'pseudospectral leapfrog', 'k-space corrected', 'Location', 'southeast'); title(sprintf('Numerical phase speed, CFL %.2f', CFL));
fprintf('\nFD2 at 10 PPW: %.3f %% (approx -(k dx)^2 (1 - CFL^2)/24 = %.3f %%)\n', 100*(2*asin(CFL*sin(pi/10))/(CFL*2*pi/10) - 1), -100*(2*pi/10)^2*(1 - CFL^2)/24);
What to look for: the measured arrival-time error of each scheme matches the dispersion-relation prediction; at 8 PPW and CFL 0.3 FD2 is late by about 2.3 % of the travel time (more than two periods over 100 λ), FD4 early by about 0.06 %, and the k-space scheme is exact to rounding.
21.14Self-check
Q1What is the value of h on axis of a disc while the arrival circle is inside the aperture, and what pressure waveform results?
\(h = c\) (\(\Theta = 2\pi\)). The pressure is \(\rho_0 c\, [v(t - z/c) - v(t - \sqrt{z^2 + a^2}/c)]\): the plane wave plus an inverted edge wave.
Q2At what depth does a 4 × 6 mm element at 3 MHz enter its far field on axis by the "SIR shorter than a period" criterion?
\((a^2 + b^2)/(2zc) = 1/f_0\) gives \(z = (a^2 + b^2) f_0/(2c) = 13 \times 3/(2 \times 1.54)\) = 12.7 mm; the point-source formula is good (≥ 0.9) only beyond about 3 times that.
Q3How much phase error does FD2 accumulate over 100 wavelengths at 10 PPW, CFL 0.3?
The phase speed is 1.5 % low, so 0.015 × 100 = 1.5 cycles (540°). FD4 at the same grid accumulates under 0.1 cycle.
Q4A two-way echo spectrum has −6 dB points at 3.4 and 6.6 MHz. Give \(f_c\), FBW and the expected −6 dB envelope length if the pulse is Gaussian.
\(f_c\) = 5.0 MHz, \(\mathrm{FBW} = 3.2/5\) = 64 %; \(\tau_{-6} = 0.883/3.2\ \mathrm{MHz}\) = 0.28 µs.
Q5What peak loss does a 0.2 mm diameter hydrophone cause at a 7.5 MHz, F# 1.5 focus?
\(\lambda F\# = 0.205 \times 1.5\) = 0.308 mm, \(b\) = 0.1 mm, \(b/\lambda F\# = 0.325\); loss \(\approx (\pi \times 0.325)^2/16 = 0.065\), about −0.58 dB (exact disc average ≈ −0.57 dB).
Q6A 3 MHz beam has \(p_r\) = 2.4 MPa in water at 4 cm. What is the derated MI at that depth?
Derating 0.3 × 3 × 4 = 3.6 dB → factor 0.661, \(p_{\mathrm{r.3}}\) = 1.59 MPa, \(\mathrm{MI} = 1.59/\sqrt{3} = 0.92\).
Q7Which simulator would you pick for (a) a 10⁶-scatterer B-mode phantom, (b) a skull aberration study, (c) the lateral modes of a 1-3 composite element?
(a) Field II (linear, fast). (b) k-Wave or another full-wave code with the heterogeneous medium (elastic if shear conversion matters). (c) Explicit FE such as PZFlex/OnScale.
Q8Why does crosstalk hurt a phased array more than a linear array?
It multiplies the element pattern by \(1 + 2\kappa \cos(kp \sin\theta)\), which is largest at broadside and smallest near \(kp \sin\theta = \pi\). Phased arrays steer to ±45° and beyond, where that factor costs several dB at −20 dB crosstalk; linear arrays mostly fire near broadside.
21.15Figures from the literature
Published figures in which the tools described above were actually run, compared with each other or with hydrophone measurements. Links only; no figures are reproduced here.
k-Wave accuracy versus points per wavelength, and k-Wave against FOCUS
Look at: Fig. 2 (transmission and reflection coefficients at a 10 % step in sound speed and density versus grid points per wavelength, for the k-space method and for first- and fourth-order finite differences: under 1 % transmission error needs about 3 PPW for k-space, 6 for FD4 and 14 for the first-order scheme), Fig. 1 (phase error after 50 wavelengths versus the reference speed in the k-space operator, zero at \(c_{\mathrm{ref}}\) = \(c_0\)) and Fig. 5 (beam of a 6 mm × 6 mm piston at 3 MHz in an absorbing medium, k-Wave at 4 PPW against the FOCUS fast near-field method). Why: the published counterpart of Demo 21.3 and of the 6 mm, 3 MHz element used throughout 21.2–21.3; note that Fig. 5 is code-to-code validation, step (3) of the ladder in 21.5.
B. E. Treeby, J. Jaros, A. P. Rendell and B. T. Cox, "Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method," J. Acoust. Soc. Am. 131(6), 4324–4336, 2012. doi:10.1121/1.4712021
k-Wave against a membrane-hydrophone scan of a clinical linear array
Look at: Fig. 1 (the tank: L9-4/38 probe, 0.4 mm PVDF membrane hydrophone, computer-controlled positioning, the layout of Figure 21.B), Fig. 2 (measured and simulated axial waveforms and spectra at 10, 20 and 30 mm in water and in olive oil) and Fig. 3 (fundamental and second-harmonic beam patterns at 0° and 20° steering). Why: step (4) of the validation ladder. The measured plane nearest the probe was used as the simulation source, which removes the transducer model from the comparison, and the authors had to apply an extra factor because the hydrophone sensitivity appeared about three times higher in oil than in water: the instrument has to be modelled before the comparison means anything.
K. Wang, E. Teoh, J. Jaros and B. E. Treeby, "Modelling nonlinear ultrasound propagation in absorbing media using the k-Wave toolbox: experimental validation," Proc. IEEE Int. Ultrasonics Symp., 523–526, 2012. doi:10.1109/ULTSYM.2012.0130
FOCUS time-space decomposition: error maps against a converged reference
Look at: Fig. 2 (transient field of a 0.6 mm-radius piston driven at 2.5 MHz, computed with a 100,000-abscissa reference integral), Fig. 3 (normalized error of the time-space decomposition with only 6 Gauss abscissas, peak 3.85 %) and Figs. 4–5 (error and computation time versus abscissas, and against the other methods compared, including the impulse response approach). Why: the near-face accuracy problem of 21.2–21.3, quantified: how many integration points a transient near field needs for a stated error.
J. F. Kelly and R. J. McGough, "A time-space decomposition method for calculating the nearfield pressure generated by a pulsed circular piston," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 53(6), 1150–1159, 2006. doi:10.1109/TUFFC.2006.1642513 (free author manuscript: PMC2117898)
Needle hydrophones: measured directivity and effective size versus frequency
Look at: Figs. 1–2 (measured directivity of 200 µm and 500 µm needle hydrophones with fits for the effective radius; compare with the rigid-disc card in 21.9, drawn for the same 0.2 and 0.5 mm diameters) and Figs. 4–5 (normalized difference between effective and geometric radius versus \(ka_g\), data from several laboratories against rigid-piston and unbaffled models). Why: at low frequency the effective diameter can exceed the geometric one by more than 100 % (502 µm on average for a 75 µm needle at 1 MHz), with deviations from rigid-piston theory below \(ka_g\) ≈ 0.7. The spatial-averaging correction of Demo 21.6 must use the effective size, not the catalogue size.
K. A. Wear, C. R. Fury and A. V. Alvarenga, "Spatial resolution limits for needle hydrophones from 0.5 to 20 MHz with implications for transcranial ultrasound," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 72(11), 1489–1496, 2025. doi:10.1109/TUFFC.2025.3610361 (NPL eprint 10365)
Field II against a hydrophone: an apodized transmit beam
Look at: Fig. 9 (five-level PWM transmit apodization: excitation sequences, ideal Gaussian and rectangular weightings, the lateral beam simulated in Field II, and the lateral beam measured with a hydrophone) and Fig. 8 (a 4.8 MHz tone burst: drive sequence, simulated and hydrophone-measured pressure). Why: a linear simulator compared with a measured field for the same array and drive, the comparison 21.8 recommends before fitting model parameters to scans.
P. R. Smith, D. M. J. Cowell and S. Freear, "Width-modulated square-wave pulses for ultrasound applications," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 60(11), 2244–2256, 2013. doi:10.1109/TUFFC.2013.6644730
21.16Go deeper
- J. A. Jensen and N. B. Svendsen, "Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 39, 262–267 (1992).
- J. A. Jensen, "Field: A program for simulating ultrasound systems," Med. Biol. Eng. Comput. 34 (Suppl. 1, Part 1), 351–353 (1996). The Field II user guide describes the sub-element and sampling approximations.
- J. A. Jensen, "A model for the propagation and scattering of ultrasound in tissue," J. Acoust. Soc. Am. 89, 182–190 (1991).
- P. R. Stepanishen, "Transient radiation from pistons in an infinite planar baffle," J. Acoust. Soc. Am. 49, 1629–1638 (1971).
- R. J. McGough, "Rapid calculations of time-harmonic nearfield pressures produced by rectangular pistons," J. Acoust. Soc. Am. 115, 1934–1941 (2004); J. F. Kelly and R. J. McGough, "A time-space decomposition method for calculating the nearfield pressure generated by a pulsed circular piston," IEEE TUFFC 53, 1150–1159 (2006). The FOCUS papers.
- B. E. Treeby and B. T. Cox, "k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields," J. Biomed. Opt. 15, 021314 (2010); B. E. Treeby, J. Jaros, A. P. Rendell and B. T. Cox, "Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method," J. Acoust. Soc. Am. 131, 4324–4336 (2012).
- M. Tabei, T. D. Mast and R. C. Waag, "A k-space method for coupled first-order acoustic propagation equations," J. Acoust. Soc. Am. 111, 53–63 (2002).
- G. R. Harris, "Progress in medical ultrasound exposimetry," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52, 717–736 (2005).
- K. A. Wear, C. R. Fury and A. V. Alvarenga, "Spatial resolution limits for needle hydrophones from 0.5 to 20 MHz with implications for transcranial ultrasound," IEEE TUFFC 72(11), 1489–1496 (2025).
- B. E. Treeby and B. T. Cox, "Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian," J. Acoust. Soc. Am. 127, 2741–2748 (2010).
- IEC 62127-1, -2, -3; IEC 62359; IEC 60601-2-37; IEC 61391-1, -2; FDA, "Marketing Clearance of Diagnostic Ultrasound Systems and Transducers" (2019, revised 2023). Titles as given in 21.7 and 21.8.