Part II · Module 10

Matched filtering & pulse compression

Coded excitation as an SNR budget tool. The matched filter makes output SNR depend on energy alone, so a long, low-amplitude code can carry the energy of a pulse the mechanical index would forbid and still compress to the resolution of its bandwidth. This module derives that result with its factor-2 conventions, prices every sidelobe fix in dB, shows what attenuation, motion and Doppler do to chirps and Golay pairs, and ends where a 510(k) reviewer starts: at fixed MI, coding raises ISPTA and TI by exactly the compression gain.

You will be able to

Derive the matched filter from Cauchy–Schwarz and state its output SNR in the real (\(2E/N_0\)) and complex (\(E/N_0\)) conventions
Predict the −6 dB width, peak sidelobe and SNR cost of a tapered chirp from \(T\), \(B\) and the window
Use Barker and Golay codes, size a least-squares mismatched filter, and compute the Golay residual for a given motion between firings
Explain range–Doppler coupling of an LFM and why it is a micron-scale bias in ultrasound
Quantify what frequency-dependent attenuation does to compression, and choose fixed, depth-adapted or equalizing filters
Convert a compression gain into penetration depth, and into MI, ISPTA.3 and TI consequences for an acoustic-output submission

Assumes

Module 03 · Convolution, LTI systems and filters: convolution, deconvolution and the Wiener filter
Module 04 · Sampling, Fourier & LTI systems: Fourier pairs, windows, sinc sidelobes, the time–bandwidth product
Module 05 · Matrices as operators: the convolution matrix and the adjoint
Module 07 · Random vectors, covariance & correlation: Gaussian noise, SNR; Module 08 · Analytic signal & IQ demodulation: analytic signal, IQ, attenuation downshift
Optional further reading: Module 11 (correlation-based delay estimation is matched filtering against a reference), Module 12 (detection and estimation bounds), Module 22 (MI, ISPTA, TI and derating, the output limits that motivate coded excitation)

Exact formulaclosed form plotted as isSimulationnumerical run in this page (FFT compression of sampled codes, noise Monte Carlo, least-squares filter solve, numerical ambiguity integral) or in the Octave scripts belowApprox. modelsimplified physics (narrowband Doppler, phase-only motion, linear-in-f attenuation without dispersion, linear acoustics); assumptions under the badgeSchematicdrawing for intuition, not computedNo transducer, Field II or acoustic-output measurement is shown: the transducer is a Gaussian band-pass, and acoustic-output numbers are linear-acoustics calculations from the stated peak pressure.

10.1The matched filter: SNR depends only on energy

An echo \(s(t - t_0)\) of known shape arrives in additive white Gaussian noise \(n(t)\). We pass the record \(x(t)\) through a linear filter \(h(t)\) and look at the output at \(t_0\). Which \(h\) maximizes the output SNR? Write the signal output and noise power in the frequency domain:

\[ y_s(t_0) = \int H(f)\,S(f)\,e^{j2\pi f t_0}\,df, \qquad \sigma_y^{2} = \frac{N_0}{2}\int \lvert H(f)\rvert^{2}\,df \]

(10.1)

\(y_s(t_0)\): filter output sampled at the echo arrival time \(t_0\), V
\(\sigma_y^{2}\): output noise variance, V²
\(S(f)\): spectrum of the echo shape \(s(t)\)
\(H(f)\): filter frequency response
\(f\): frequency, MHz
\(N_0/2\): two-sided power spectral density of the real white noise, W/Hz (or V²/Hz)
\(t_0\): echo arrival time, s
\(j\): imaginary unit

\[ \mathrm{SNR} = \frac{\left\lvert \int H S\, e^{j2\pi f t_0}\,df \right\rvert^{2}}{\frac{N_0}{2}\int \lvert H\rvert^{2}\,df} \;\le\; \frac{\int \lvert S\rvert^{2}\,df}{N_0/2} = \frac{2E}{N_0} \]

(10.2)

\(\mathrm{SNR}\): peak output signal power over output noise variance, linear
\(H,\ S\): filter response and echo spectrum, as in (10.1)
\(f,\ t_0,\ j\): frequency, arrival time and imaginary unit, as in (10.1)
\(E\): echo energy \(\int s^2\,dt = \int\lvert S\rvert^2 df\) (Parseval), J or V²·s
\(N_0\): one-sided noise power spectral density; the two-sided level is \(N_0/2\)

Cauchy–Schwarz: \(\lvert\langle u, v\rangle\rvert^2 \le \lVert u\rVert^2\lVert v\rVert^2\) with \(u = H\), \(v = S^* e^{-j2\pi f t_0}\).

\[ \text{equality} \iff H(f) = \kappa\, S^*(f)\, e^{-j2\pi f t_0} \iff h(t) = \kappa\, s^*(t_0 - t) \]

(10.3)

\(\kappa\): arbitrary nonzero gain (does not change SNR)
\(H(f),\ S(f)\): filter response and echo spectrum, as in (10.1)
\(S^*\): complex conjugate of \(S\)
\(h(t)\): matched-filter impulse response
\(s(t)\): echo shape; \(s^*\) is its complex conjugate (equal to \(s\) for real RF)
\(t_0,\ j\): arrival time and imaginary unit, as in (10.1)
\(\iff\): if and only if

The matched filter is the time-reversed, conjugated echo. Its output is the cross-correlation of the record with the template: \(y(t) = \int x(u)\, s^*(u - t + t_0)\,du\), the correlation receiver.

Two consequences drive the rest of the module. First, the maximum SNR contains only \(E\): the shape of \(s\) is free. A 0.3 µs pulse and a 20 µs chirp with the same energy give the same peak SNR, so duration can buy energy while bandwidth buys resolution. Second, the optimum is a ratio of inner products, so it is flat near its peak: a moderately wrong filter costs little SNR. For a Gaussian echo spectrum of std \(\sigma_s\) (amplitude) and a Gaussian filter of std \(\sigma_h\) offset by \(\Delta f\),

\[ \eta = \frac{\left\lvert\int H S\right\rvert^{2}}{\int\lvert H\rvert^{2}\int\lvert S\rvert^{2}} = \frac{2\sigma_s\sigma_h}{\sigma_s^{2} + \sigma_h^{2}}\, \exp\!\left(-\frac{\Delta f^{2}}{\sigma_s^{2} + \sigma_h^{2}}\right) \]

(10.4)

\(\eta\): SNR efficiency relative to the matched filter, \(0 < \eta \le 1\)
\(H,\ S\): filter response and echo spectrum, as in (10.1)
\(\sigma_s,\ \sigma_h\): standard deviations of the echo and filter amplitude spectra, MHz
\(\Delta f\): centre-frequency offset between filter and echo, MHz

With \(r = \sigma_h/\sigma_s\): \(2r/(1 + r^2)\). \(r\) = 2 or ½ costs 0.97 dB; \(\Delta f = \sigma_s\) with \(r\) = 1 costs 2.17 dB. Narrowband approximation (it ignores the negative-frequency lobe).

Colored noise. If the noise PSD is \(N(f)\) rather than flat, the same argument with \(u = H\sqrt{N}\) gives \(H = S^* e^{-j2\pi f t_0}/N(f)\): whiten, then match. The maximum becomes \(\int \lvert S(f)\rvert^2/N(f)\,df\). This is also the reason a filter matched to the electrical drive waveform is wrong: the template must be the received echo, after both passes through the transducer.

Data and noise model	Signal energy	Noise	Peak SNR of the matched filter
Real RF, continuous	\(E = \int s^2\,dt\)	two-sided PSD \(N_0/2\)	\(2E/N_0\)
Real RF, sampled at \(f_s\)	\(\sum s_k^2 = E f_s\)	\(\sigma^2 = (N_0/2) f_s\) per sample	\(\sum s_k^2/\sigma^2 = 2E/N_0\)
Analytic signal or IQ, \(\mathrm{SNR} = \lvert y\rvert^2/\mathrm{E}\lvert n_{\mathrm{out}}\rvert^2\)	\(\int\lvert s_a\rvert^2 = 2E\) (scale-free ratio)	both quadratures counted	\(E/N_0\) (3.01 dB lower)
Analytic signal, SNR of \(\mathrm{Re}\{y\}\) at the peak	as above	in-phase quadrature only	\(2E/N_0\) again
Complex white noise, discrete, \(\mathrm{E}\lvert n_k\rvert^2 = \sigma^2\)	\(\sum\lvert s_k\rvert^2\)	\(\sigma^2\) total (\(\sigma^2/2\) per quadrature)	\(\sum\lvert s_k\rvert^2/\sigma^2\)
Communications \(E_b/N_0\)	energy per bit	one-sided PSD \(N_0\)	quoted as \(E_b/N_0\); the peak SNR is \(2E_b/N_0\)

Demo 10.1

Matched vs mismatched filtering of a pulse in noise

SimulationSampled Gaussian-spectrum pulse in white Gaussian noise, filtered in JS; Monte Carlo over 4000 noise records. The dashed level is the closed form √(2E/N₀)·√η.

record (signal + noise)echo / noise-free outputfilter output with noise√(2E/N₀) · √η

Input 2E/N₀ = Σs²/σ² 10 dB

Filter / echo bandwidth r 1.00

Filter centre offset Δf 0.0 MHz

Try: at \(r\) = 1 and \(\Delta f\) = 0 the filter is matched and the output peak reaches \(\sqrt{2E/N_0}\) noise standard deviations. Halve or double the bandwidth (\(r\) = 0.5 or 2): the cost is about 1 dB. Push \(\Delta f\) to 2 MHz (\(1.6\sigma\)) and the loss exceeds 5 dB. At 10 dB the echo barely stands out of the record (its best single sample is only 3.5 dB above the noise, because the energy is spread over many samples), yet the output peak sits at 3.2 noise standard deviations.

Ultrasound link

You already use matched filters. DAS is the spatial matched filter of a point target (Module 05: the adjoint, not the inverse). Correlation-based delay estimation is matched filtering against a reference A-line (Module 11), and its peak is the maximum-likelihood delay in white noise. The receive band-pass filter on any scanner approximates a matched filter to the transducer's two-way response; the result above says that getting it within a factor 2 in bandwidth is good enough for SNR. TGC multiplies signal and noise alike and changes no SNR.

Factor-2 traps

Most "3 dB discrepancies" between a simulation and \(2E/N_0\) are convention, not physics: complex noise counted in both quadratures, one-sided versus two-sided \(N_0\), or energy of the analytic signal (\(2E\)) used with the noise of the real one. State the convention next to every SNR you report, and measure SNR as mean²/variance of the filter output at the peak sample over many noise realizations.

10.2Time–bandwidth product and compression gain

A conventional pulse couples energy to resolution. Its energy is about \(p_{\mathrm{peak}}^2\tau/(\rho c)\) per unit area and its axial resolution is \(c\tau/2\) with \(\tau \approx 1/B\). The mechanical index caps \(p_{\mathrm{peak}}\), so with a short pulse the only way to more energy is a longer pulse, which destroys resolution. Pulse compression breaks the coupling: transmit a code of duration \(T \gg 1/B\) with bandwidth \(B\), then compress it with the matched filter to a mainlobe of width \(\approx 1/B\).

\[ \text{compression ratio} = \frac{T}{1/B} = TB, \qquad G = \frac{E_{\mathrm{code}}}{E_{\mathrm{ref}}} \ \text{ at equal peak amplitude} \]

(10.5)

\(T\): code duration, µs
\(B\): swept or occupied bandwidth, MHz
\(TB\): time-bandwidth product, dimensionless
\(G\): compression (SNR) gain, linear; quote in dB as \(10\log_{10}G\)
\(E_{\mathrm{code}},\ E_{\mathrm{ref}}\): energy of the code and of the reference short pulse, J

The matched-filter SNR is proportional to energy, (10.2), so the SNR gain is the energy ratio at the same peak pressure.

\[ \begin{gathered} \text{rect reference of length } 1/B:\ G = TB; \\ \text{Gaussian reference, same −6 dB band:}\ G = \frac{T}{\sqrt{\pi}\,\sigma_t} = \frac{2\sqrt{\pi}}{2.3548}\,TB = 1.505\,TB \end{gathered} \]

(10.6)

\(G\): compression gain, as in (10.5)
\(T,\ B\): code duration (µs) and bandwidth (MHz), as in (10.5)
\(\sigma_t\): Gaussian envelope std, \(\sigma_t = 1/(2\pi\sigma_f)\), µs
\(\sigma_f\): amplitude-spectrum std, \(\sigma_f = B/2.3548\) for a −6 dB width \(B\), MHz
\(2.3548\): \(2\sqrt{2\ln 2}\), ratio of −6 dB amplitude width to std

Constant-envelope code of amplitude \(A\): energy \(A^2T\) (complex envelope). Gaussian envelope \(A\,e^{-t^2/2\sigma_t^2}\): energy \(\sqrt\pi\,\sigma_t A^2\). TB = 50 gives 18.77 dB; a Monte Carlo with \(10^4\) trials gives 18.70 dB (m10_chirp_compression.m).

The mechanism is phase alignment. A chirp delays each frequency by a different group delay, \(\tau_g(f) = (f - f_0)T/B\) for an up-chirp. The matched filter has the opposite group delay, so after filtering all frequencies arrive at the same instant: the output spectrum \(\lvert S(f)\rvert^2\) is real and zero-phase, and its transform is as short as the bandwidth allows. A phase code does the same thing with a random-looking phase spectrum that the filter conjugates away.

Intuition

Think of the code as a deliberately defocused pulse in time, exactly like a defocused aperture in space. The matched filter applies the conjugate delay law, the way a focusing law cancels the quadratic phase across the aperture (Module 04). An LFM is literally a quadratic phase in time, and compression is focusing in range.

10.3The linear-FM chirp: mainlobe, sidelobes, tapers

The LFM chirp sweeps its instantaneous frequency linearly through the band:

\[ s(t) = w\!\left(\frac{t}{T}\right)\cos\!\left(2\pi\left(f_0 t + \frac{B}{2T}t^{2}\right)\right), \quad \lvert t\rvert \le \frac{T}{2}, \qquad f_i(t) = f_0 + \frac{B}{T}\,t \]

(10.7)

\(s(t)\): transmitted chirp, normalized amplitude
\(t\): time from the centre of the chirp, µs
\(w\): amplitude envelope of argument \(t/T\) (rect for an unweighted chirp)
\(T,\ B\): chirp duration (µs) and swept bandwidth (MHz)
\(f_0\): centre frequency, MHz
\(f_i(t)\): instantaneous frequency, MHz
\(B/T\): sweep rate, MHz/µs

In complex baseband: \(\tilde s(t) = w(t/T)\,e^{j\pi(B/T)t^2}\).

\[ \lvert S(f)\rvert \approx \sqrt{\frac{T}{B}}\; w\!\left(\frac{(f - f_0)T}{B}\right), \qquad \angle S(f) \approx -\pi\frac{(f - f_0)^{2}\,T}{B} + \frac{\pi}{4} \]

(10.8)

\(S(f)\): chirp spectrum (stationary-phase approximation, valid for \(TB \gg 1\))
\(\angle S\): spectral phase, rad
\(f,\ f_0\): frequency and centre frequency, MHz
\(T,\ B\): chirp duration (µs) and bandwidth (MHz), as in (10.7)
\(w\): envelope shape of (10.7), here mapped onto frequency

Stationary phase, valid for \(TB \gg 1\): each frequency is emitted at the instant \(t = (f - f_0)T/B\), so the time envelope maps onto the spectrum. For finite \(TB\) the band edges carry Fresnel-integral ripple.

\[ \frac{\lvert R(\tau)\rvert}{R(0)} = \left(1 - \frac{\lvert\tau\rvert}{T}\right)\left\lvert \mathrm{sinc}\!\left(B\tau\left(1 - \frac{\lvert\tau\rvert}{T}\right)\right)\right\rvert, \quad \lvert\tau\rvert < T; \qquad w_{-6\,\mathrm{dB}} \to \frac{1.207}{B},\ \ \mathrm{PSL} \to -13.26\ \mathrm{dB} \]

(10.9)

\(R(\tau)\): matched-filter output (autocorrelation) at lag \(\tau\); \(R(0)\) is its peak
\(\tau\): lag, µs
\(T,\ B\): chirp duration (µs) and bandwidth (MHz)
\(\mathrm{sinc}(x)\): \(\sin(\pi x)/(\pi x)\)
\(w_{-6\,\mathrm{dB}}\): full width at half amplitude of the compressed envelope, µs
\(\mathrm{PSL}\): peak sidelobe level relative to the mainlobe, dB (20 log)

Exact for the rect-envelope complex chirp (m10_chirp_compression.m matches it to \(3\times10^{-5}\)). As \(TB \to \infty\) it becomes \(\mathrm{sinc}(B\tau)\): −6 dB width \(1.207/B\), first null \(1/B\), peak sidelobe −13.26 dB. At TB = 10: \(1.234/B\) and −14.9 dB; at TB = 200: \(1.208/B\) and −13.3 dB.

With \(B\) = 5 MHz the compressed −6 dB width is 241 ns, an axial resolution of \(c \cdot 1.207/(2B)\) = 0.186 mm, independent of \(T\). The −13 dB range sidelobes extend over \(\pm T\), i.e. \(\pm cT/2\) in depth (±7.7 mm for 10 µs), and are unacceptable in a B-mode image displayed over 50–60 dB. They are lowered by amplitude tapering, and each placement has a price, all relative to the rect chirp at the same peak amplitude:

\[ \text{receive taper (mismatched):}\ \eta = \frac{\left(\sum w\right)^{2}}{N\sum w^{2}}; \qquad \text{transmit taper:}\ \eta = \overline{w^{2}}; \qquad \text{split } \sqrt{w} \text{ on both sides:}\ \eta = \overline{w} \]

(10.10)

\(\eta\): SNR efficiency (loss factor) relative to the untapered code at equal peak amplitude
\(w\): taper samples across the code, peak 1
\(N\): number of samples in the code
\(\sum\): sum over the \(N\) samples
\(\overline{\,\cdot\,}\): mean over the code

Hann: −1.76, −4.26 and −3.01 dB. The transmit loss is an energy loss at fixed peak pressure; if MI is not the binding limit you can raise the amplitude and recover it.

Window, TB = 50	Receive: PSL / width·B / loss	Transmit	Split √w	Ideal spectral window PSL
rect	−13.5 dB / 1.21 / 0.00 dB	same	same	−13.3 dB
Tukey 25 %	−13.6 / 1.37 / −0.42	−13.7 / 1.43 / −0.74	−13.7 / 1.37 / −0.58	n/a
Hann	−31.3 / 2.00 / −1.76	−46.3 / 2.59 / −4.26	−31.1 / 2.00 / −3.01	−31.5
Hamming	−36.5 / 1.82 / −1.34	−49.9 / 2.41 / −4.01	−38.9 / 1.82 / −2.68	−42.7
Blackman	−35.8 / 2.30 / −2.37	−81.4 / 3.11 / −5.16	−53.3 / 2.30 / −3.77	−58.1

Read the Hamming and Blackman receive rows: they miss their textbook levels by 6 and 22 dB. The compressed output is the transform of the window times the chirp spectrum, and the spectrum is only approximately a rect: its Fresnel ripple has sidelobes of its own that the window does not touch. A transmit taper smooths the envelope, and with it the band edges where the ripple lives, which is why its PSL goes much lower. Nonlinear FM (a sweep rate that slows near the band centre) shapes the spectrum by dwell time instead of amplitude and avoids the mismatch loss; it is the radar answer and has been studied for ultrasound as well.

Demo 10.2

LFM chirp: duration, bandwidth, taper and where you put it

SimulationSampled chirp compressed by FFT correlation in JS; widths and PSL measured on the result.

transmitted RF and compressed enveloperect chirp, matched (reference)

Duration T 10 µs

Swept bandwidth B 5 MHz

Taper

Applied on

Try: set the taper to rect and sweep \(T\) at fixed \(B\): the mainlobe stays at \(1.2/B\) while the sidelobes stay near −13 dB and stretch over \(\pm T\). Then compare Hamming on receive with Hamming on transmit: transmit reaches about −50 dB but costs 4 dB of energy at the same peak pressure. Reduce \(TB\) below 20 and watch the Fresnel ripple raise every tapered sidelobe floor.

10.4Phase codes: Barker, mismatched filters and Golay pairs

A binary phase code transmits \(N\) chips of duration \(T_c\) with phases 0 or \(\pi\) (amplitudes \(c_k = \pm 1\)). In ultrasound a chip is typically one or a few carrier cycles. The matched filter output, sampled at chip spacing, is the aperiodic autocorrelation

\[ R(k) = \sum_i c_i\, c_{i+k}, \qquad R(0) = N, \qquad G = N \ \text{ at equal peak amplitude} \]

(10.11)

\(R(k)\): autocorrelation of the code at lag \(k\)
\(c_i\): chip amplitudes, ±1
\(i,\ k\): chip index and lag in chips
\(N\): code length in chips
\(G\): compression gain, as in (10.5)

The mainlobe is one chip wide (resolution set by \(T_c\), i.e. by bandwidth); the gain is the energy ratio \(N\).

Barker codes have \(\lvert R(k)\rvert \le 1\) for every \(k \ne 0\), the smallest possible for a binary code. Barker-13 (+ + + + + − − + + − + − +) has sidelobes of 0 or +1, so \(\mathrm{PSL} = 20\log_{10}(1/13) = -22.28\) dB and \(G = 10\log_{10}13 = 11.14\) dB. No Barker code of odd length above 13 exists and no longer one is known, so binary codes need help. A least-squares mismatched filter \(\mathbf{h}\) of length \(L\) minimizes \(\sum_{k\ne0}\lvert(c * h)_k\rvert^2\) for a unit peak:

\[ \mathbf{h} = \arg\min_{\mathbf{h}} \lVert \mathbf{C}\mathbf{h} - \mathbf{d}\rVert^{2} = (\mathbf{C}^{T}\mathbf{C})^{-1}\mathbf{C}^{T}\mathbf{d}, \qquad \text{loss} = \frac{\lvert (c * h)_0\rvert^{2}}{\lVert\mathbf{h}\rVert^{2}\,\lVert\mathbf{c}\rVert^{2}} \]

(10.12)

\(\mathbf{h}\): mismatched filter coefficients, length \(L\)
\(\mathbf{C}\): \((N+L-1)\times L\) convolution matrix of the code (Module 05)
\(\mathbf{c}\): code vector, \(N\) chips
\(\mathbf{d}\): centred spike (desired output)
\(N,\ L\): code length and filter length in chips
\(*\): linear convolution; \((c * h)_0\) is the output at the peak lag
\(\lVert\cdot\rVert\): Euclidean norm
\((\cdot)^{T}\): transpose
\(\arg\min\): argument that minimizes the cost
\(\text{loss}\): SNR loss of the mismatched filter relative to the matched one, linear (≤ 1)

Barker-13: \(L\) = 13 → −24.0 dB, \(L\) = 39 → −38.5 dB, \(L\) = 65 → −58.6 dB, for 0.14–0.21 dB of SNR loss.

Golay complementary pairs \((A, B)\) have autocorrelations whose sidelobes cancel exactly:

\[ R_A(k) + R_B(k) = 2N\,\delta(k), \qquad (A, B) \to (A\,|\,B,\ A\,|\,{-B}) \ \text{ doubles } N \]

(10.13)

\(A,\ B\): the two complementary codes of a Golay pair
\(R_A(k),\ R_B(k)\): autocorrelations of \(A\) and \(B\) at lag \(k\) (chips)
\(N\): length of each code in chips
\(\delta(k)\): Kronecker delta, 1 at \(k = 0\) and 0 elsewhere
\(A\,|\,B\): concatenation of \(A\) and \(B\), giving a pair of length \(2N\)

Fire A, compress with A; fire B, compress with B; add. Gain \(2N\) over the two firings (\(N\) = 16: 15.05 dB), sidelobes exactly zero for a static target, also after any band-limiting that is identical for both firings (the RF simulation leaves −102 dB, numerical floor). \(A\,|\,B\) is concatenation.

The price is two firings and a static target. If the target moves axially by \(d\) between firings, echo B arrives with an extra round-trip phase \(\varphi = 4\pi f_0 d/c\) (for \(d\) much smaller than a chip). Then \(R_A(k) + e^{j\varphi}R_B(k) = R_A(k)(1 - e^{j\varphi})\) off the peak and \(2N\cos(\varphi/2)\,e^{j\varphi/2}\) at the peak:

\[ \mathrm{PSL}_{\mathrm{residual}} = \frac{\max_{k\ne0}\lvert R_A(k)\rvert}{N}\,\tan\frac{\varphi}{2}, \qquad \varphi = \frac{4\pi f_0 d}{c} \]

(10.14)

\(\mathrm{PSL}_{\mathrm{residual}}\): residual peak sidelobe level relative to the peak, linear
\(R_A(k),\ N\): autocorrelation and code length, as in (10.13)
\(k\): lag in chips
\(d\): axial displacement between the two firings, µm
\(f_0\): centre frequency, MHz
\(c\): sound speed, 1540 m/s
\(\varphi\): round-trip phase shift of echo B, rad

\(N\) = 16 (\(\max\lvert R_A\rvert\) = 5) at 5 MHz: \(d\) = 2 µm → −37.9 dB, 5 µm → −29.9 dB, 20 µm → −17.4 dB. At a 5 kHz PRF, 2 µm per firing is 10 mm/s of axial velocity. The 1-cycle-chip RF simulation agrees within 0.5 dB.

SchematicTransmit chain drawn for orientation; each lane's waveform is computed from its formula but the three lanes use different time scales (10 µs, 2.6 µs, two firings one PRI apart), so they are not to a common scale.

Coded excitation from code memory to compression, and three transmit sequences. Top: an LFM sweeps 2.5 → 7.5 MHz in 10 µs (TB = 50), so echoes from the first \(cT/2\) = 7.7 mm overlap the transmission. Middle: Barker-13 with one 5 MHz cycle per chip, phase flipped by π on the − chips (2.6 µs, gain 11.1 dB). Bottom: a Golay pair is two firings; the receiver adds the two compressed outputs, which cancels the sidelobes only if the target has not moved between them, (10.13)–(10.14).

Demo 10.3

Barker-13, least-squares mismatched filter, Golay pairs and motion

SimulationApprox. modelChip-level correlations and least-squares filter solve in JS; motion enters as a phase-only rotation of echo B (valid for d much smaller than a chip).

Code and filter

Golay length N

LS filter length L 39

Axial motion between firings d 2.0 µm

Try: with the Golay pair at \(d\) = 0 every sidelobe vanishes (drawn at the −80 dB floor). Move \(d\) to 5 µm (25 mm/s at 5 kHz PRF): the residual is the single-code sidelobe pattern scaled by \(\tan(\varphi/2)\). For Barker, slide \(L\) from 13 to 79 and watch the PSL fall while the SNR loss stays near 0.2 dB.

Ultrasound link

Golay pairs suit slow targets. Combined with multi-pulse contrast schemes (pulse inversion, amplitude modulation) they multiply the number of firings, and every term that differs between firings, motion or a nonlinear bubble response, appears as residual sidelobes. For plane-wave and ULM acquisitions the second firing halves the effective frame rate. A chirp compresses in one firing.

10.5The ambiguity function: Doppler tolerance and range–Doppler coupling

A scatterer moving at radial velocity \(v\) compresses or stretches the echo. In the narrowband approximation this is a Doppler shift \(\nu = 2vf_0/c\), and the matched filter output for a shifted echo is a cut through the ambiguity function

\[ \chi(\tau, \nu) = \int \tilde s(t)\, \tilde s^*(t - \tau)\, e^{j2\pi\nu t}\,dt, \qquad \iint \lvert\chi(\tau, \nu)\rvert^{2}\,d\tau\,d\nu = E^{2} \]

(10.15)

\(\chi(\tau,\nu)\): ambiguity function of the code
\(\tilde s\): complex envelope of the transmitted code (\(\tilde s^*\) its conjugate)
\(t\): time, s
\(\tau\): delay mismatch, s
\(\nu\): Doppler shift \(2vf_0/c\), Hz
\(E\): energy of \(\tilde s\)
\(j\): imaginary unit

Conventions differ in sign and in where the conjugate sits. The volume property says ambiguity can be moved around the \((\tau, \nu)\) plane, never removed.

\[ \text{LFM:}\quad \frac{\lvert\chi(\tau, \nu)\rvert}{E} = \left(1 - \frac{\lvert\tau\rvert}{T}\right)\left\lvert \mathrm{sinc}\!\left(\left(\nu + \frac{B\tau}{T}\right)(T - \lvert\tau\rvert)\right)\right\rvert \]

(10.16)

\(\chi(\tau,\nu),\ E\): ambiguity function and code energy, as in (10.15)
\(\tau,\ \nu\): delay mismatch (µs) and Doppler shift (MHz), as in (10.15)
\(T,\ B\): chirp duration (µs) and bandwidth (MHz)
\(\mathrm{sinc}\): \(\sin(\pi x)/(\pi x)\), as in (10.9)

The ridge runs along \(\nu = -B\tau/T\). A Doppler shift \(\nu\) moves the compressed peak by \(\lvert\Delta\tau\rvert \approx \lvert\nu\rvert T/B\) and lowers it by about \(20\log_{10}(1 - \lvert\nu\rvert/B)\). At finite \(TB\) the \((1 - \lvert\tau\rvert/T)\) factor pulls the peak back toward zero by \(3/(\pi^2B^2T)\) (1.2 ns for \(T\) = 10 µs, \(B\) = 5 MHz; checked against a numerical maximization). The sign flips between up- and down-chirps.

In radar the coupling matters because \(\nu\) can approach \(B\). In ultrasound \(v/c \approx 10^{-3}\): at 1 m/s and 5 MHz, \(\nu\) = 6.49 kHz. For \(T\) = 10 µs and \(B\) = 5 MHz the ridge position \(\nu T/B\) is 13.0 ns and the actual peak, pulled back by the \((1 - \lvert\tau\rvert/T)\) envelope, sits at 11.8 ns: 9 µm of apparent depth, about 5 % of the 241 ns mainlobe, for a loss of 0.01 dB. The SNR loss is irrelevant; the shift is a velocity-dependent bias of the kind that matters for speckle tracking and microbubble localization, and it cancels when up- and down-chirps are averaged. The narrowband model itself holds while \((2v/c)\,TB \ll 1\) (0.065 in this example). Phase codes have a thumbtack-like ambiguity near \(\nu = 0\) whose Doppler width is \(\approx 1/T\); since \(\nu T\) = 0.065 here, they lose essentially nothing either.

Demo 10.4

Ambiguity maps and the ultrasound Doppler band

SimulationApprox. modelAmbiguity integral evaluated numerically in JS on sampled envelopes; narrowband (Doppler-shift) model of motion.

±ν_D = ±2vf₀/c (f₀ = 5 MHz)|χ| in dB, −30 to 0

Waveform (T = 10 µs)

Scatterer velocity v 1.0 m/s

Try: on the LFM map the ultrasound band is a hairline around \(\nu = 0\) even at 5 m/s, because the map has to span \(\pm B\) to show the ridge. Switch to the unmodulated burst: its Doppler resolution is \(1/T\), but its range response is the full 10 µs triangle, which is why nobody images with it. Barker-13 concentrates the volume in a narrow spike plus a sidelobe floor.

10.6What tissue does to compression: attenuation, band limits, nonlinearity

The received echo from depth \(z\) is the transmitted code filtered twice by the transducer and once by the round-trip attenuation. For attenuation linear in frequency, \(\alpha\) in dB/cm/MHz and \(\alpha_0 = \alpha/8.686\) in Np/cm/MHz,

\[ E(f; z) = S(f)\, H_{\mathrm{tr}}(f)\, e^{-2\alpha_0 f z} \]

(10.17)

\(E(f;z)\): amplitude spectrum of the echo from depth \(z\)
\(S(f)\): spectrum of the transmitted pulse
\(H_{\mathrm{tr}}(f)\): two-way transducer response
\(\alpha_0\): attenuation slope, Np/cm/MHz (\(\alpha\) in dB/cm/MHz ÷ 8.686)
\(f,\ z\): frequency in MHz, depth in cm (path \(2z\))

The power spectrum carries \(e^{-4\alpha_0 f z}\). Module 08 derives the Gaussian-spectrum downshift \(-2\alpha_0\sigma_f^2 z\) (amplitude std \(\sigma_f\)); a chirp's flat band shifts less.

Compress with the filter matched to the surface echo, \(H = E^*(f; 0)\). The output spectrum \(E(f; z)E^*(f; 0) = \lvert S\rvert^2\lvert H_{\mathrm{tr}}\rvert^2 e^{-2\alpha_0 f z}\) is real and non-negative, so in this zero-phase model the compressed envelope still peaks at lag 0. What attenuation takes is bandwidth: the mainlobe widens, the spectral taper becomes asymmetric so sidelobes rise, and Cauchy–Schwarz charges a mismatch loss

\[ \eta(z) = \frac{\left(\int \lvert E_0\rvert\,\lvert E_z\rvert\,df\right)^{2}}{\int \lvert E_0\rvert^{2}\,df \int \lvert E_z\rvert^{2}\,df} \]

(10.18)

\(\eta(z)\): SNR efficiency of the fixed (surface-matched) filter at depth \(z\), \(0 < \eta \le 1\)
\(E_0,\ E_z\): echo spectra \(E(f;0)\) and \(E(f;z)\) of (10.17)
\(f\): frequency, MHz

\(E_0 = E(f; 0)\), \(E_z = E(f; z)\). 2.5–7.5 MHz chirp, \(T\) = 10 µs, 80 % Gaussian transducer, 0.5 dB/cm/MHz: centroid 5.00 → 3.68 MHz at 6 cm; loss 0.26, 0.97, 2.01, 3.23, 4.55 dB at 2, 4, 6, 8, 10 cm; −6 dB width 338 → 374 ns; PSL −29.8 → −26.1 dB.

The uncompressed echo behaves differently: an LFM maps frequency to time, so weighting toward low frequencies moves its energy centroid by about \((f_c - f_0)T/B\), −2.5 µs at 6 cm. Do not read depth from the raw chirp envelope. Three receive strategies (filters written as transfer functions applied to the echo):

Fixed (matched at \(z\) = 0), \(H = E_0^*\): simplest, loses SNR and sidelobe level with depth.
Depth-adapted matched \(H = E_z^*\): optimal SNR by definition, but its output spectrum \(\lvert E_z\rvert^2\) is narrower still, so the mainlobe widens (475 ns at 6 cm, −20 dB width 1.5 µs).
Equalizing (Wiener-type) \(H = E_z^*\lvert E_0\rvert^2/(\lvert E_z\rvert^2 + \epsilon)\): restores the surface pulse shape (376 ns at 6 cm) and pays 5.8 dB relative to the adapted filter, because it re-amplifies the attenuated high frequencies and their noise.

\[ H_W(f) = \frac{S^*(f)}{\lvert S(f)\rvert^{2} + N(f)/P(f)} \]

(10.19)

\(H_W(f)\): Wiener-type filter response
\(S^*(f),\ S(f)\): conjugate of the echo spectrum and the spectrum itself
\(N(f)\): noise power spectrum
\(P(f)\): power spectrum of the reflectivity (signal prior)

The Wiener compression filter interpolates between the matched filter (\(N/P\) large: maximum SNR, wide mainlobe) and the inverse filter (\(N/P \to 0\): shortest pulse, noise amplified). Least-squares mismatched filters (10.12) are its sidelobe-targeted cousins.

Transducer bandwidth acts as a spectral taper and, through the LFM time–frequency map, as a transmit taper: with an 80 % two-way Gaussian response the surface PSL is already −29.8 dB without any explicit window. Part of the swept band is also lost, so the realized gain is below \(10\log_{10}(1.505\,TB)\) computed from the electrical sweep. Nonlinear propagation at diagnostic pressures moves energy to harmonics that the fundamental filter treats as loss. For quadratic nonlinearity the second harmonic of \(a(t)\cos(\omega t + \phi(t))\) carries \(2\phi(t)\): a chirp's harmonic is a chirp with twice the centre frequency and twice the sweep, compressible with its own matched filter, while a biphase code's phases 0 and \(\pi\) both become 0 mod \(2\pi\), so its harmonic is an uncoded long pulse. Near field: echoes from depths below \(cT/2\) overlap the transmission and are not fully compressed (7.7 mm for 10 µs), so systems shorten the code near the probe.

Demo 10.5

Chirp through attenuating tissue: fixed, adapted and equalizing filters

SimulationApprox. modelFFT compression in JS of a sampled chirp × Gaussian two-way transducer × linear-in-f attenuation without dispersion (same model as m10_attenuation_mismatch.m).

|E(f; 0)||E(f; z)| / fixed filter outputdepth-adaptedequalizer

Attenuation α 0.50 dB/cm/MHz

Depth z 6.0 cm

Transducer −6 dB band (two-way) 80 %

Try: at \(\alpha\) = 0 all three filters coincide. Raise \(z\): the fixed filter's mainlobe widens and its loss grows roughly quadratically at first; the adapted filter's loss stays at 0 by definition while its −20 dB width explodes; the equalizer holds the width and pays in SNR. The compressed peak never moves; the raw echo centroid does.

Pitfall

A PSL number alone can flatter a filter. The depth-adapted filter at 6 cm shows −35.8 dB PSL because its broad mainlobe skirt has no null for the sidelobe search to start from. Always report the −20 dB (or −40 dB) width or an integrated sidelobe level next to the PSL.

10.7The SNR budget, penetration and acoustic output

In a noise-limited region the echo SNR falls by the two-way attenuation, \(2\alpha f_0 z\) dB (narrowband). A compression gain \(G\) dB moves the noise-limited depth by

\[ \Delta z = \frac{G}{2\,\alpha\, f_0} \quad [\mathrm{cm}] \]

(10.20)

\(\Delta z\): added noise-limited depth, cm
\(G\): compression gain, dB
\(\alpha\): attenuation slope, dB/cm/MHz (one way)
\(f_0\): centre frequency, MHz

15 dB at 0.5 dB/cm/MHz and 5 MHz: 3.0 cm. TB-50 chirp (18.77 dB): 3.75 cm at 5 MHz, 6.26 cm at 3 MHz, 2.50 cm at 7.5 MHz. Barker-13 (11.14 dB): 2.23 cm at 5 MHz.

The rule is conservative for broadband pulses. With an adapted filter the SNR at depth is \(\int\lvert E_0(f)\rvert^2 e^{-4\alpha_0 f z}df\), and the deep echo is dominated by the low end of the band where the loss per cm is smaller. For a 60 % Gaussian spectrum at 5 MHz and 0.5 dB/cm/MHz, truncated at its −24 dB points (2 and 8 MHz), a 40 dB short-pulse budget reaches 9.7 cm instead of 8.0 cm, and 15 dB of coding adds 4.9 cm instead of 3.0 cm, with a centre frequency of 2.8 MHz at the new depth; without the truncation the Gaussian tail keeps feeding low frequencies and the numbers become 9.8 cm, 5.7 cm and 2.1 MHz (m10_snr_budget_output.m).

Equal energy, equal SNR. In white noise, \(K\) averaged short pulses and one code with \(K\) times their energy give the same SNR: one TB-50 chirp matches 75 averaged pulses (Monte Carlo 14.36 vs 14.12 dB). Coding does not create energy; it delivers it in one firing at a peak pressure the short pulse could not use.

Acoustic output. That last sentence is the regulatory content of coded excitation. The indices respond to different properties of the waveform:

\[ \begin{gathered} \mathrm{MI} = \frac{p_{r.3}}{\sqrt{f_{\mathrm{awf}}}}\ \ (\mathrm{MPa},\ \mathrm{MHz}), \qquad \mathrm{PII} = \int \frac{p^{2}(t)}{\rho c}\,dt, \\ I_{\mathrm{SPTA}} = \mathrm{PII} \times (\text{pulses per second at the point}), \qquad \mathrm{TI} \propto \text{acoustic power} \end{gathered} \]

(10.21)

\(\mathrm{MI}\): mechanical index, dimensionless
\(p_{r.3}\): peak rarefactional pressure derated by 0.3 dB/cm/MHz, at the point where the derated pulse-intensity integral is largest, MPa
\(f_{\mathrm{awf}}\): acoustic working frequency, MHz
\(\mathrm{PII}\): pulse-intensity integral, J/m²
\(p(t)\): instantaneous pressure waveform, Pa
\(\rho,\ c\): density (kg/m³) and sound speed (m/s); \(\rho c \approx 1.5\) MRayl
\(I_{\mathrm{SPTA}}\): spatial-peak temporal-average intensity, mW/cm²
\(\mathrm{TI}\): thermal index, dimensionless

FDA Track 3 global maxima (derated): \(I_{\mathrm{SPTA.3}} \le\) 720 mW/cm², and MI ≤ 1.9 or \(I_{\mathrm{SPPA.3}} \le\) 190 W/cm²; MI and TI are displayed to the user.

At the same peak pressure and frequency, MI is unchanged while PII, \(I_{\mathrm{SPTA}}\) and power rise by the energy ratio, which is the compression gain. Worked example (m10_snr_budget_output.m): MI = 1.9 at 5 MHz means \(p_{r.3}\) = 4.25 MPa; a 60 % Gaussian pulse at 200 pulses/s at a point gives \(I_{\mathrm{SPTA.3}}\) = 27 mW/cm²; a TB-50 chirp at the same peak gives 2003 mW/cm² (×75.2, 18.8 dB), well above 720. Either the pulse rate drops to 72 Hz or the peak falls by ×0.60 (MI 1.14). The available headroom at that point was 720/27 = 26.7, i.e. 14.3 dB of coding gain before \(I_{\mathrm{SPTA}}\) binds. Beyond that, the energy limit binds, and at a binding energy limit coding gains nothing against white noise. The constant-amplitude chirp's \(I_{\mathrm{SPPA.3}} \approx p^2/(2\rho c)\) is about 600 W/cm², so it complies through the MI branch of the Track 3 rule, not the \(I_{\mathrm{SPPA}}\) branch.

Demo 10.6

SNR budget and acoustic output of a coded mode

Approx. modelNarrowband rule and broadband Gaussian-spectrum integral for depth; linear acoustics, plane-wave intensity p²/(ρc), no nonlinear propagation, for the output indices.

narrowband ruleshort pulse (broadband)coded (broadband)Track 3 I_SPTA.3 limit

Centre frequency f₀ 5.0 MHz

Attenuation α 0.50 dB/cm/MHz

Chirp TB (B = 0.6 f₀) 50

Receive taper

Short-pulse SNR at z = 0 40 dB

MI (both waveforms) 1.90

Pulses per second at the point 200 Hz

Try: at the defaults the coded mode reaches 3.8 cm deeper (narrowband rule) but exceeds the \(I_{\mathrm{SPTA.3}}\) limit by 2.8×. Lower the MI until the coded bar clears 720 and note how much gain is left. Raise the pulse rate toward ultrafast values (several kHz at the same point): even the short pulse approaches the limit, and coding must run at low MI.

510(k) angle

For a new coded mode, the acoustic-output section is measured in that mode at its worst-case settings (longest code, highest pulse rate, maximum drive, every focal setting), not inferred from the short-pulse mode. MI may be unchanged at the same \(p_r\), but \(I_{\mathrm{SPTA.3}}\) and the thermal indices scale with energy, so expect TI to become the displayed index that moves. Two details need a stated method: the acoustic working frequency of a swept waveform (the zero-crossing and −3 dB spectral definitions can disagree for a chirp, and MI scales as \(1/\sqrt{f_{\mathrm{awf}}}\)), and nonlinear propagation at high drive, which changes both \(p_r\) and the spectrum in water. The applicable standards are IEC 60601-2-37 and IEC 62359 (index determination), with the FDA diagnostic-ultrasound guidance defining the Track 3 limits.

10.8Visual summary

Thirteen questions an expert answers in five seconds, each with the picture computed from its formula; the badge says whether it is a closed form, a numerical run or a simplified model.

10.9How experts think about it

Experts treat coded excitation as accounting. The matched filter turns SNR into energy, so every design question becomes: where does the energy come from (duration, at a peak pressure capped by MI), what does it cost (\(I_{\mathrm{SPTA}}\), TI, frame rate, blind zone), and how much leaks into range sidelobes (taper, mismatch, Fresnel ripple, attenuation, motion)? They quote gains in dB against a stated reference, losses against the matched filter, and depth gains as \(G/(2\alpha f_0)\), and they remember that coding fights noise, not clutter.

Peak SNR depends only on the echo energy: 2E/N₀, whatever the waveform.

\(\mathrm{SNR}_{\max} = 2E/N_0\)\(h(t) = s^*(t_0 - t)\)

Cauchy–Schwarz on \(\lvert\int HS\rvert^2\) against the noise power \((N_0/2)\int\lvert H\rvert^2\) gives the bound \(\int\lvert S\rvert^2/(N_0/2) = 2E/N_0\), reached only by \(H \propto S^*\), (10.2)–(10.3). Shape does not appear, so duration (energy) and bandwidth (resolution) can be chosen separately: that freedom is pulse compression. m10_matched_filter_snr.m: a 5 MHz pulse at \(\sum s^2/\sigma^2\) = 10.00 dB measures 9.94 dB over 8000 trials.

SimulationSampled pulse and chirp of equal energy, discrete correlation in JS.

White Gaussian noise, known waveform, output read at the true delay. Colored noise needs prewhitening (\(H = S^*/N(f)\)); clutter and reverberation are not noise.

State your SNR convention: real RF gives 2E/N₀, the complex |y|² definition gives E/N₀.

factor 2 = 3.01 dB\(\sigma^2 = (N_0/2)\,f_s\)\(E_b/N_0\): one-sided \(N_0\)

Sampled real RF with two-sided PSD \(N_0/2\): \(\sigma^2 = N_0 f_s/2\) and \(\sum s^2 = E f_s\), so \(\sum s^2/\sigma^2 = 2E/N_0\). The analytic signal has noise in both quadratures; \(\lvert y\rvert^2/\mathrm{E}\lvert n\rvert^2\) counts both and reads \(E/N_0\), while \(\mathrm{Re}\{y\}\) at the peak reads \(2E/N_0\). The script prints 9.94, 6.95 and 9.94 dB for the three definitions on the same data.

SimulationBars are the 8000-trial Monte Carlo results of m10_matched_filter_snr.m (Octave).

This is bookkeeping, not physics. The real cost of the quadrature noise appears only in envelope detection at low SNR.

Matched filters are flat near the optimum: a factor-2 bandwidth error costs about 1 dB.

\(\eta = 2r/(1 + r^2)\)r = 2 or ½: −0.97 dBoffset 1σ: −2.17 dB

For Gaussian echo and filter spectra with stds \(\sigma_s\), \(\sigma_h\) and offset \(\Delta f\), the efficiency is \(\frac{2\sigma_s\sigma_h}{\sigma_s^2+\sigma_h^2}\exp\!\left(-\frac{\Delta f^2}{\sigma_s^2+\sigma_h^2}\right)\), (10.4). The first factor is symmetric in \(r\) and \(1/r\) and has zero slope at \(r\) = 1. This is why one receive filter serves many depths for SNR purposes, and why attenuation mismatch shows up first as sidelobes, not as lost SNR.

Approx. modelNarrowband Gaussian spectra, negative-frequency lobe ignored.

Narrowband Gaussian spectra. For broadband filters (\(r\) = 3) the exact loss is −2.55 dB against −2.22 dB predicted. For coded waveforms the sidelobe consequences of mismatch matter more than the SNR.

Compression gain is the time–bandwidth product, measured at equal peak amplitude.

\(G \approx TB\)vs Gaussian, same −6 dB band: \(1.505\,TB\)TB 50: 18.8 dB

MI caps the peak pressure, and SNR follows energy. A constant-envelope code of amplitude \(A\) has energy \(A^2T\); a Gaussian pulse \(\sqrt\pi\,\sigma_t A^2\) with \(\sigma_t = 2.3548/(2\pi B)\). The ratio is \(2\sqrt\pi\,TB/2.3548 = 1.505\,TB\), (10.6); against a rect pulse of length \(1/B\) it is exactly \(TB\). Monte Carlo at TB = 50: 18.70 dB against 18.77 dB predicted.

Exact formulaThe single point is the Octave Monte Carlo result.

White noise only. Tapers subtract their loss, the transducer removes part of the swept band, and attenuation mismatch subtracts more at depth. Always say which reference pulse the gain is quoted against.

Compressed resolution is set by B, not by T: the −6 dB width is 1.207/B.

B = 5 MHz: 241 ns = 0.186 mmfirst null \(\approx 1/B\)

For \(TB \gg 1\) the chirp spectrum is close to a rect of width \(B\), so the autocorrelation tends to \(T\,\mathrm{sinc}(B\tau)\), which falls to 0.5 at \(B\tau\) = 0.6034: full width \(1.207/B\). The exact rect-chirp form (10.9) gives \(1.234/B\) at TB = 10 and \(1.208/B\) at TB = 200. Axial −6 dB resolution is \(c \cdot 1.207/(2B)\).

Exact formula

Rect spectrum. Hann on receive widens it to \(2.00/B\); the transducer and attenuation shrink the effective \(B\).

An unweighted chirp has −13 dB range sidelobes spread over ±cT/2, which no B-mode image tolerates.

PSL → −13.26 dBT = 10 µs: ±7.7 mmdisplay: 50–60 dB

It is the uniform-aperture sinc of Module 04, now along depth. A specular echo 40 dB above speckle would paint its sidelobes 27 dB above the speckle over ±7.7 mm, filling any nearby cyst or vessel. Sidelobes have to sit below the displayed dynamic range relative to the brightest reflector that can appear next to a dark region, so practical targets are −40 to −50 dB or better.

Exact formula

The transducer itself tapers the spectrum: with an 80 % two-way band the surface PSL is −29.8 dB without any explicit window. Integrated sidelobe level, not only PSL, decides cyst fill-in.

A receive taper costs (Σw)²/(NΣw²); a transmit taper costs mean(w²) at fixed peak pressure.

Hann rx: −1.76 dB, PSL −31.3 dBHann tx: −4.26 dBsplit √w: −3.01 dB

A receive taper makes the filter mismatched; Cauchy–Schwarz gives the efficiency \(\lvert\sum w\rvert^2/(N\sum w^2) = 0.25/0.375 = 0.667\) for Hann, (10.10). A transmit taper keeps the filter matched but removes energy when the peak is capped: \(\overline{w^2}\) = 0.375. Split tapering pays \(\overline{w}\) = 0.5. Numbers are for TB = 50 in m10_chirp_compression.m.

Exact formula

Losses are relative to the rect chirp at the same peak. If \(I_{\mathrm{SPTA}}\) or TI binds instead of MI, raise the amplitude to restore the energy and the transmit-taper loss disappears.

Fresnel ripple puts a floor under any receive taper: Hamming gives −36.5 dB at TB 50, not −42.7.

Blackman rx: −35.8 dB (ideal −58.1)Blackman tx: −81.4 dB

The chirp spectrum is a rect only in the stationary-phase limit, (10.8). At finite \(TB\) its band edges ripple (Fresnel integrals), and the compressed output is the transform of window × rippled spectrum, so the ripple's own sidelobes leak above the window's ideal level. Tapering on transmit smooths the envelope and with it the band edges where the ripple lives, which is why the transmit rows go far lower. Larger \(TB\) and nonlinear FM reduce the problem.

SimulationChirp spectrum computed by direct numerical Fourier sum in JS (TB = 50).

Numbers for TB = 50 and the taper applied across the whole code; transducer band-limiting also smooths the edges in practice.

Barker-13 is as good as binary codes get: 13 times the energy, sidelobes at 1/13.

G = 11.14 dBPSL −22.28 dBLS filter L = 39: −38.5 dB for 0.20 dB

Barker sequences have \(\lvert R(k)\rvert \le 1\) off the peak, the minimum for ±1 chips. No Barker code of odd length above 13 exists and none longer is known. A least-squares mismatched filter \((\mathbf{C}^{T}\mathbf{C})^{-1}\mathbf{C}^{T}\mathbf{d}\) of three times the code length lowers the PSL to −38.5 dB for 0.20 dB of SNR loss, (10.12); \(L\) = 65 reaches −58.6 dB for 0.21 dB.

Exact formula

Chip-level model. One-cycle chips at \(f_0\) are heavily filtered by the transducer, and the filter must be designed on the received waveform; the second harmonic of a biphase code carries no code.

Golay pairs cancel exactly, but only if nothing moves between the two firings.

\(R_A + R_B = 2N\delta\)N 16, d = 2 µm: −37.9 dBd = 20 µm: −17.4 dB

A displacement \(d\) between firings rotates echo B by \(\varphi = 4\pi f_0 d/c\). Off the peak \(R_A(1 - e^{j\varphi})\) survives with magnitude \(2\lvert R_A\rvert\sin(\varphi/2)\); the peak becomes \(2N\cos(\varphi/2)\). The residual PSL is \((\max\lvert R_A\rvert/N)\tan(\varphi/2)\), (10.14). At a 5 kHz PRF and 5 MHz, 10 mm/s gives 2 µm and −37.9 dB for \(N\) = 16; 10 cm/s (cardiac wall) gives −17.4 dB.

Approx. modelPhase-only motion model, d much smaller than a chip.

Phase-only model, valid for \(d\) much smaller than a chip; the 1-cycle-chip RF simulation agrees within 0.5 dB. Nonlinear propagation or bubble responses that differ between firings break it too. Frame rate halves.

A chirp turns Doppler into a range shift Δτ ≈ νT/B, and in ultrasound that is microns.

1 m/s at 5 MHz: ν = 6.49 kHzT 10 µs, B 5 MHz: 12 ns ≈ 9 µmloss 0.01 dB

The LFM ambiguity ridge runs along \(\nu = -B\tau/T\), so a Doppler-shifted echo compresses near \(\tau = -\nu T/B\) (11.8 ns instead of 13.0 ns here, the finite-\(TB\) pull \(3/(\pi^2B^2T)\)) with a peak lowered by about \(20\log_{10}(1 - \nu/B)\), (10.16). With \(2v/c \approx 10^{-3}\) the loss is negligible, but the shift is about 5 % of the mainlobe: a velocity-dependent position bias for speckle tracking or bubble localization. Up- and down-chirps shift in opposite directions, so averaging them cancels it.

Approx. modelNarrowband ambiguity function of a rect LFM; peak found numerically.

Narrowband ambiguity, valid while \((2v/c)\,TB \ll 1\) (0.065 here). Only the axial velocity component produces Doppler. The sign depends on sweep direction and convention.

Attenuation does not move the compressed peak; it eats bandwidth: wider mainlobe, higher sidelobes, mismatch loss.

0.5 dB/cm/MHz, 6 cm: 5.00 → 3.68 MHzloss 2.0 dB\(w_{-6\,\mathrm{dB}}\) 338 → 374 ns

In a dispersion-free model \(e^{-2\alpha_0 f z}\) is a zero-phase filter, so \(E_zE_0^*\) is real and non-negative and the envelope peaks at lag 0 (simulated shift 0.00 ns at every depth), (10.17)–(10.18). The raw echo's centroid moves by \(\approx (f_c - f_0)T/B\), −2.5 µs at 6 cm, because an LFM maps frequency to time. The adapted filter recovers the 2 dB but widens to 475 ns; the equalizer keeps 376 ns and costs 5.8 dB.

Approx. modelIdealized 2.5–7.5 MHz band, 80 % Gaussian transducer, linear-in-f attenuation, no dispersion.

Linear-in-f, dispersion-free attenuation. Real tissue adds weak dispersion (Kramers–Kronig) and frequency-dependent scattering; Doppler, not attenuation, is what moves an LFM peak.

Coding buys G/(2αf₀) of depth: 15 dB at 0.5 dB/cm/MHz and 5 MHz is 3 cm.

\(\Delta z = G/(2\alpha f_0)\)TB 50 at 3 MHz: 6.3 cmbroadband: 4.9 cm, not 3.0

The two-way loss is \(2\alpha f_0 z\) dB with \(\alpha\) in dB/cm/MHz, so \(G\) dB of SNR moves the noise-limited depth by \(G/(2\alpha f_0)\), (10.20). For a broadband pulse with an adapted filter the deep echo is dominated by low frequencies whose loss per cm is smaller: for a 60 % Gaussian spectrum at 5 MHz cut at 2 and 8 MHz, 15 dB adds 4.9 cm, and the spectrum there is centred at 2.8 MHz (5.7 cm and 2.1 MHz if the Gaussian tail is not cut).

Approx. modelNarrowband attenuation rule, noise-limited region.

Noise-limited region, white noise, linear-in-f attenuation; the coded spectrum is flatter than a Gaussian, so the broadband correction is approximate. Clutter-limited depths do not move.

Coding helps against noise, never against clutter: reverberation and off-axis echoes scale with the transmit.

clutter-to-signal: unchangednoise-limited depth: \(+G/(2\alpha f_0)\)

Clutter, reverberation and sidelobe echoes are filtered copies of the same transmitted code, so the matched filter compresses them with the same gain as the target, and range sidelobes add an axial clutter term of their own. Only receiver noise is independent of the transmit, and only it is beaten by more energy. So coding extends the noise-limited depth (transcranial, deep abdomen, deep ULM) and does nothing for a near-field cyst limited by reverberation.

SchematicTwo bookkeeping lines (+G dB and 0 dB), not a simulation.

Where noise and clutter are comparable, lowering the noise floor still improves contrast; the gain is then smaller than \(G\).

At fixed MI, coding raises ISPTA and TI by exactly the compression gain.

MI 1.9 at 5 MHz: \(p_{r.3}\) = 4.25 MPaTB 50: ISPTA.3 ×75 (18.8 dB)Track 3: ≤ 720 mW/cm²

MI uses only the peak rarefactional pressure and the frequency, both unchanged. \(\mathrm{PII} = \int p^2/(\rho c)\,dt\) scales with energy, \(I_{\mathrm{SPTA}}\) = PII × pulse rate, and TI scales with power, (10.21). At 200 pulses/s at a point: 27 → 2003 mW/cm². The headroom 720/27 = 26.7 means up to 14.3 dB of coding is free in \(I_{\mathrm{SPTA}}\) terms there; beyond it the energy limit binds and coding gains nothing in white noise.

Approx. modelLinear acoustics, plane-wave intensity, 27 mW/cm² short-pulse value from m10_snr_budget_output.m.

Linear propagation and a nominal \(f_{\mathrm{awf}}\). Submissions use measured, derated values in the actual mode; nonlinear propagation changes \(p_r\) and the spectrum.

Equal energy, equal SNR: one TB-50 chirp matches 75 averaged short pulses.

Monte Carlo: 14.36 vs 14.12 dB1 firing vs 75

Coherently summing \(K\) matched-filtered pulses multiplies the signal amplitude by \(K\) and the noise variance by \(K\), so the SNR grows by \(K\): exactly what one code with \(K\) times the energy gives. Coding therefore buys time (frame rate, no motion during the average) and lower peak pressure, not energy. That is the right way to explain it to a regulator or a reviewer.

Exact formula

White noise independent across firings and a static target during the average; averaging has no range sidelobes, coding does.

The blind zone is cT/2: a 10 µs code hides the first 7.7 mm.

5 µs: 3.85 mm20 µs: 15.4 mm40 µs: 30.8 mm

When the same elements transmit and receive, echoes from depths below \(cT/2\) return while the code is still being transmitted, so they are blocked or only partially compressed. Systems therefore shorten the code near the probe and lengthen it with depth, which also matches where the extra SNR is needed.

Exact formula

Separate transmit and receive apertures reduce the problem; partial compression degrades gradually rather than as a hard edge.

A biphase code vanishes at the second harmonic; a chirp survives with twice the sweep.

\(2\phi \in \{0, 2\pi\}\)chirp: \(2f_0\), \(2B\) over \(T\)

Quadratic nonlinearity produces \(a^2(t)\cos(2\omega t + 2\phi(t))/2\) at the second harmonic. For a biphase code \(\phi \in \{0, \pi\}\), so \(2\phi \equiv 0\): the harmonic is an uncoded pulse of length \(T\) with no compression. For a chirp \(2\phi(t)\) is again quadratic, a chirp centred at \(2f_0\) sweeping \(2B\), compressible with a filter matched to it. That is why harmonic coded imaging uses chirps.

Exact formula

Weak quadratic nonlinearity; the squared envelope acts as an extra taper. Microbubbles respond with memory and depend on peak pressure and duration, so contrast imaging needs its own analysis.

10.10Cheat sheet

Concept	Formula	Ultrasound use
Matched filter	\(h(t) = s^(t_0 - t)\), \(H = S^ e^{-j2\pi f t_0}\)	Pulse compression, receive filter, DAS, correlation delay estimation
Peak SNR	\(2E/N_0\) (real, two-sided \(N_0/2\)); \(E/N_0\) for complex \(\lvert y\rvert^2\)	Budget calculations; check simulations against it
Colored noise	\(H = S^*/N(f)\), \(\mathrm{SNR} = \int\lvert S\rvert^2/N\,df\)	Prewhitening when the noise spectrum is shaped
Mismatch efficiency	\(\eta = \frac{\lvert\int HS\rvert^2}{\int\lvert H\rvert^2\int\lvert S\rvert^2}\); Gaussian \(\frac{2r}{1+r^2}e^{-\Delta f^2/(\sigma_s^2+\sigma_h^2)}\)	Tolerance of receive filters, tapers, attenuation
Compression gain	\(G = TB\) (rect \(1/B\) ref.), \(1.505\,TB\) (Gaussian ref., same −6 dB band)	SNR gain at the same peak pressure
LFM chirp	\(f_i = f_0 + (B/T)t\); \(\lvert R\rvert = (1-\lvert\tau\rvert/T)\lvert\mathrm{sinc}(B\tau(1-\lvert\tau\rvert/T))\rvert\)	Single-firing coded excitation
Resolution, sidelobes	\(w_{-6\,\mathrm{dB}} = 1.207/B\); \(\Delta z = 1.207c/(2B)\); PSL −13.26 dB	Axial PSF of a compressed chirp
Taper costs	rx \((\sum w)^2/(N\sum w^2)\); tx \(\overline{w^2}\); split \(\overline{w}\)	Hann: −1.76 / −4.26 / −3.01 dB
Barker-13	\(\lvert R(k{\ne}0)\rvert \le 1\); PSL −22.28 dB; \(G\) = 11.14 dB	Short binary codes; LS mismatched filter for −40 to −60 dB
Golay pair	\(R_A + R_B = 2N\delta\); residual \((\max\lvert R_A\rvert/N)\tan(2\pi f_0 d/c)\)	Two-firing coding; motion sensitivity
LFM Doppler coupling	\(\Delta\tau \approx \nu T/B\) (\(TB \gg 1\)), \(\nu = 2vf_0/c\); valid for \((2v/c)TB \ll 1\)	Velocity-dependent position bias, cancelled by up/down chirps
Attenuated echo	\(E = S H_{\mathrm{tr}} e^{-2\alpha_0 f z}\), \(\alpha_0 = \alpha/8.686\)	Fixed, adapted or equalizing compression filters
Wiener compression	\(H = S^*/(\lvert S\rvert^2 + N/P)\)	Trade SNR for mainlobe width and sidelobes
Penetration gain	\(\Delta z = G/(2\alpha f_0)\) cm	Depth bought by coding (narrowband, noise-limited)
Acoustic output	\(\mathrm{MI} = p_{r.3}/\sqrt{f_{\mathrm{awf}}}\); \(I_{\mathrm{SPTA}}\) = PII × rate; TI ∝ power	At fixed MI, \(I_{\mathrm{SPTA}}\) and TI scale by \(G\); Track 3: \(I_{\mathrm{SPTA.3}} \le\) 720 mW/cm², MI ≤ 1.9
Blind zone	\(z_{\min} = cT/2\)	Shorter codes in the near field

10.11Interview and qualifying-exam questions

I1Derive the matched filter and its output SNR. · screen

Signal output at \(t_0\) is \(\int HS\,e^{j2\pi f t_0}df\), noise variance \((N_0/2)\int\lvert H\rvert^2 df\) for real white noise of two-sided PSD \(N_0/2\). Cauchy–Schwarz bounds the ratio by \(\int\lvert S\rvert^2/(N_0/2) = 2E/N_0\), with equality for \(H = \kappa S^* e^{-j2\pi f t_0}\), i.e. \(h(t) = \kappa s^*(t_0 - t)\). The output is the cross-correlation of the record with the template; only \(E\) appears, not the shape.

Follow-up: what if the noise is colored? Whiten and match: \(H = S^* e^{-j2\pi f t_0}/N(f)\), \(\mathrm{SNR} = \int\lvert S\rvert^2/N(f)\,df\).

I2Your simulation shows matched-filter SNR exactly 3 dB below 2E/N₀. Bug? · screen

Probably a convention. If the data are analytic or IQ and SNR is defined as \(\lvert y\rvert^2/\mathrm{E}\lvert n_{\mathrm{out}}\rvert^2\), both noise quadratures are counted and the result is \(E/N_0\). Taking \(\mathrm{Re}\{y\}\) at the peak restores \(2E/N_0\) (the script measures 6.95 vs 9.94 dB). Other usual causes: one-sided \(N_0\), or the analytic signal's energy \(2E\) paired with the real signal's noise.

Follow-up: does envelope detection really lose 3 dB? Not at high SNR, where the quadrature noise adds only a second-order term to the envelope; the loss appears as a detection-threshold effect at low SNR.

I3Why use a chirp instead of simply driving the probe harder? · onsite

Peak rarefactional pressure is capped by the MI limit (1.9, Track 3), by pulser voltage and by nonlinear distortion at high pressure. SNR depends on energy, so at a fixed peak the only lever is duration; compression restores the \(1/B\) resolution. A TB-50 chirp gives 18.8 dB against a Gaussian pulse of the same band, about 3.8 cm more depth at 5 MHz and 0.5 dB/cm/MHz.

Follow-up: what limits \(T\)? Blind zone \(cT/2\), \(I_{\mathrm{SPTA}}\)/TI (energy scales with \(T\)), frame rate in multi-line or plane-wave sequences, attenuation mismatch, and transmit/receive switching.

I4Your chirp shows −13 dB range sidelobes. How do you get to −50 dB, and what does it cost? · onsite

Taper. On receive alone, Hann gives −31.3 dB for 1.76 dB and a 1.65× wider mainlobe, but Hamming and Blackman stall near −36 dB at TB = 50 because of Fresnel ripple. On transmit the PSL drops to −46 (Hann) or −81 dB (Blackman) at a 4–5 dB energy cost at fixed peak; split tapering sits between. Alternatives: nonlinear FM (spectrum shaped by sweep rate, no mismatch loss), least-squares or Wiener mismatched filters, and larger \(TB\).

Follow-up: why does the transmit taper beat the ideal-window number for Hann? Because the output spectrum is then closer to the window squared, and the smooth envelope suppresses the band-edge ripple.

I5Golay pair or chirp for deep plane-wave ULM at 3 MHz? · research

Golay gives exactly zero static sidelobes but needs two firings per angle, halving the compounded frame rate. With \(N\) = 16 and 5 cm/s axial motion at a 5 kHz PRF, \(d\) = 10 µm, \(\varphi = 4\pi \cdot 3\ \mathrm{MHz} \cdot 10\ \mu\mathrm{m}/1540\ \mathrm{m/s}\) = 0.245 rad and the residual is −28.3 dB; tissue and bubble motion at higher speed degrade it further, and the nonlinear bubble echo differs between firings. A tapered chirp compresses in one firing, survives motion (about 12 ns of Doppler bias per m/s at 5 MHz, \(B\) = 5 MHz, \(T\) = 10 µs), and its harmonic remains compressible. For ULM I would choose a tapered chirp and budget its sidelobes against the bubble-to-tissue ratio after SVD filtering.

Follow-up: how would you rescue Golay? Estimate the inter-firing displacement (e.g. Kasai phase, Module 11) and rotate echo B before summing.

I6What does frequency-dependent attenuation do to a compressed chirp, and how do you fix it? · research

It tapers the received spectrum toward low frequencies: centre 5.00 → 3.68 MHz at 6 cm for 0.5 dB/cm/MHz with an 80 % transducer. With the surface filter the peak does not move (zero-phase model), but the mainlobe widens (338 → 374 ns), the PSL rises (−29.8 → −26.1 dB) and the mismatch loss is 2.0 dB. A depth-adapted filter recovers the 2 dB but widens to 475 ns; an equalizer keeps 376 ns and costs 5.8 dB. Depth-dependent filter banks interpolate between them.

Follow-up: why does the raw chirp envelope move by −2.5 µs while the compressed peak stays? The LFM maps frequency to time, so a spectral tilt shifts the envelope centroid by \((f_c - f_0)T/B\); the matched filter removes the quadratic phase first.

I7Is range–Doppler coupling of an LFM a concern in ultrasound? · research

Not for SNR, sometimes for position. \(\Delta\tau \approx \nu T/B\) with \(\nu = 2vf_0/c\): at 1 m/s, 5 MHz, \(T\) = 10 µs, \(B\) = 5 MHz the ridge gives 13.0 ns and the exact peak 11.8 ns (finite-\(TB\) correction \(3/(\pi^2B^2T)\)), i.e. 9 µm of apparent depth for 0.01 dB of loss. That is small against the 0.19 mm resolution but not against ULM localization precision or µm-scale elastography displacements, and it depends on velocity. Alternate up- and down-chirps, or average them, to cancel it.

Follow-up: when does the narrowband ambiguity function fail? When \((2v/c)\,TB\) approaches 1, i.e. when the echo's time compression over the code length is comparable to the resolution: never at tissue speeds with clinical codes.

I8You add a coded B-mode to a cleared system. What changes in the 510(k) acoustic-output section? · onsite

It must be measured in the new mode at worst-case settings. At the same peak pressure MI is unchanged, but PII, \(I_{\mathrm{SPTA.3}}\) and power scale by the energy ratio, so a TB-50 chirp at MI 1.9 and 200 pulses/s at a point goes from 27 to 2003 mW/cm², over the 720 mW/cm² Track 3 limit; TIs rise by the same factor. You would cap MI (×0.60 to 1.14 here) or the pulse rate for the coded mode, define \(f_{\mathrm{awf}}\) for a swept waveform, and check derated values with nonlinear propagation in water.

Follow-up: can you claim the SNR gain if \(I_{\mathrm{SPTA}}\) already binds? No: at a binding energy limit the gain against white noise is zero; coding helps where MI, not energy, was the binding constraint.

I9How much depth does 15 dB of compression gain buy? · screen

\(\Delta z = G/(2\alpha f_0)\): 3.0 cm at 0.5 dB/cm/MHz and 5 MHz, 5.0 cm at 3 MHz. For a broadband pulse with an adapted filter it is more, 4.9 cm at 5 MHz with a 60 % Gaussian band cut at its −24 dB points, because the deep spectrum downshifts to where attenuation per cm is lower.

Follow-up: and in a clutter-limited region? Nothing: clutter scales with the transmit.

I10Will coded excitation improve the contrast of a liver cyst at 6 cm in an obese patient? · onsite

Only to the extent the cyst is noise-limited. Reverberation from the abdominal wall and off-axis clutter are copies of the transmitted code and compress with the same gain, and range sidelobes of bright interfaces add axial fill-in. If the noise floor sits above the clutter at that depth, coding lowers it by up to \(G\) and contrast improves; otherwise harmonic imaging or better aberration and clutter suppression matter more.

Follow-up: how would you tell which regime you are in? Turn the transmit off (or receive with transmit muted) and compare the noise floor with the signal level inside the cyst at the same gain.

I11What is the PSL of Barker-13, and how would you improve it without a longer code? · screen

Sidelobes are 0 or 1 against a peak of 13: −22.28 dB, with an 11.14 dB gain. A least-squares mismatched filter of length 39 gives −38.5 dB for 0.20 dB of SNR loss (65 taps: −58.6 dB). Concatenating or nesting codes, or switching to Golay pairs, are the other options.

Follow-up: why not design the filter on the ±1 chip sequence directly? Because the received chips are band-limited by the transducer; the mismatched filter must be designed on the received waveform.

10.12MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run. The numbers quoted below are what the scripts print.

Download the scripts: m10_matched_filter_snr.m · m10_chirp_compression.m · m10_barker_golay.m · m10_attenuation_mismatch.m · m10_snr_budget_output.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m10_matched_filter_snr.mMatched filter SNR in three conventions, Gaussian bandwidth and centre mismatch · change bw, EsN0_dB, r_list, Ntrial · runtime ≈ 6 s

%% m10_matched_filter_snr.m
% Matched filter output SNR in white Gaussian noise, three conventions, and
% the cost of using the wrong filter.
%   - Real RF samples, white noise of variance sigma^2: peak SNR = sum(s.^2)/sigma^2,
%     which is 2E/N0 for a two-sided noise PSD N0/2 (sigma^2 = (N0/2)*fs).
%   - Complex (analytic / IQ) data, SNR defined as |y|^2 / E|n_out|^2: 3 dB lower (E/N0).
%   - Gaussian band-pass filter of the wrong bandwidth (ratio r) or centre:
%     efficiency 2r/(1+r^2) and exp(-D^2/(sf^2+sh^2)) (narrowband theory).
% What to change: bw (fractional bandwidth), EsN0_dB, r_list, Ntrial.
% Expected output: measured RF SNR within 0.1 dB of theory; analytic |y|^2 SNR
% 3.0 dB lower; Re{} of the analytic output back at the RF value; r = 0.5 costs 0.97 dB.
% Runtime about 6 s in Octave.
clear; close all; rng(1);

%% Pulse and noise
fs = 40e6;              % sampling rate (Hz)
f0 = 5e6;               % centre frequency (Hz)
bw = 0.60;              % -6 dB fractional bandwidth of the amplitude spectrum
sf = bw*f0/2.3548;      % std of the Gaussian amplitude spectrum (Hz)
st = 1/(2*pi*sf);       % std of the Gaussian envelope (s)
t  = (-3e-6:1/fs:3e-6)';
s0 = exp(-t.^2/(2*st^2)).*cos(2*pi*f0*t);
EsN0_dB = 10;                            % target sum(s.^2)/sigma^2 in dB
sigma = 1;
s = s0*sqrt(10^(EsN0_dB/10)*sigma^2/sum(s0.^2));
SNRth = sum(s.^2)/sigma^2;               % = 2E/N0
fprintf('Pulse: %d samples, sum(s^2)/sigma^2 = %.2f dB (theory 2E/N0)\n', numel(s), 10*log10(SNRth));
fprintf('Gain over the best single sample: %.2f dB\n', 10*log10(sum(s.^2)/max(s.^2)));

%% Monte Carlo: matched filter on RF and on the analytic signal (and mismatched filters)
Nchunk = 2; Nper = 4000; Ntrial = Nchunk*Nper;  % 8000 trials in chunks to save memory
Nrec = 512; i0 = 256;                          % record length, echo centre index
x0 = zeros(Nrec,1); L = numel(s); x0(i0-(L-1)/2 : i0+(L-1)/2) = s;
win = i0-(L-1)/2 : i0+(L-1)/2;
wA = zeros(Nrec,1); wA(1) = 1; wA(2:Nrec/2) = 2; wA(Nrec/2+1) = 1;   % analytic-signal weights
sA = ifft(fft(x0).*wA); sA = sA(win);
r_list = [0.25 0.5 0.75 1 1.5 2 3 4];
effTh = 2*r_list./(1+r_list.^2); effEx = zeros(size(r_list));
Hs = cell(size(r_list)); rows = zeros(size(r_list)); Nf = Nrec + 512;
for k = 1:numel(r_list)
    sth = st/r_list(k);                          % filter envelope std
    th = (-ceil(6*sth*fs):ceil(6*sth*fs))'/fs;
    h = exp(-th.^2/(2*sth^2)).*cos(2*pi*f0*th);  % symmetric, so h(-t) = h(t)
    ys = conv(s, h); [~, kp] = max(abs(ys));
    effEx(k) = ys(kp)^2/(sigma^2*sum(h.^2))/SNRth;   % exact deterministic efficiency
    Hs{k} = fft(h, Nf); rows(k) = win(1) + kp - 1;  % output index of the signal peak
end
yRF = zeros(1,Ntrial); yA = zeros(1,Ntrial); yK = zeros(numel(r_list),Ntrial);
for c = 1:Nchunk
    X = repmat(x0,1,Nper) + sigma*randn(Nrec,Nper);
    idx = (c-1)*Nper + (1:Nper);
    yRF(idx) = s' * X(win,:);                    % matched filter output at the peak
    XA = ifft(fft(X).*repmat(wA,1,Nper));
    yA(idx) = sA' * XA(win,:);                   % sA' is the conjugate transpose
    Xf = fft([X; zeros(Nf-Nrec,Nper)]);
    for k = 1:numel(r_list)
        Y = real(ifft(Xf.*repmat(Hs{k},1,Nper)));
        yK(k,idx) = Y(rows(k),:);
    end
end
SNRrf = mean(yRF)^2/var(yRF);
mA = mean(yA);
SNRc  = abs(mA)^2/mean(abs(yA-mA).^2);
SNRre = real(mA)^2/var(real(yA));
effMC = (mean(yK,2).^2./var(yK,0,2))'/SNRth;
fprintf('\nMatched filter, %d trials\n', Ntrial);
fprintf('  RF (real), coherent        : %.2f dB  (theory %.2f)\n', 10*log10(SNRrf), 10*log10(SNRth));
fprintf('  analytic, |y|^2/E|n|^2     : %.2f dB  (theory %.2f = E/N0)\n', 10*log10(SNRc), 10*log10(SNRth/2));
fprintf('  analytic, Re{y} (coherent) : %.2f dB\n', 10*log10(SNRre));

fprintf('\nGaussian filter of bandwidth ratio r (filter/echo):\n   r    theory 2r/(1+r^2)   exact     Monte Carlo\n');
for k = 1:numel(r_list)
    fprintf('%5.2f   %7.2f dB        %7.2f dB  %7.2f dB\n', r_list(k), 10*log10(effTh(k)), 10*log10(effEx(k)), 10*log10(effMC(k)));
end
for D = [0.5 1 2]
    h = exp(-t.^2/(2*st^2)).*cos(2*pi*(f0+D*sf)*t);
    ys = conv(s,h); Nz = 2048; wz = [1; 2*ones(Nz/2-1,1); 1; zeros(Nz/2-1,1)];
    A = abs(ifft(fft(ys, Nz).*wz));               % envelope of the filter output
    fprintf('centre offset %.1f sf (%.2f MHz): loss %.2f dB, theory exp(-D^2/2) = %.2f dB\n', D, D*sf/1e6, ...
        10*log10(max(A)^2/(sigma^2*sum(h.^2))/SNRth), 10*log10(exp(-D^2/2)));
end

%% Figures
figure('Position',[100 100 900 560]);
subplot(2,1,1);
tt = (0:Nrec-1)'/fs*1e6; X = repmat(x0,1,1) + sigma*randn(Nrec,1);
plot(tt, X(:,1), 'Color', [0.6 0.6 0.6]); hold on; plot(tt, x0, 'b', 'LineWidth', 1.5);
xlim([2 11]); xlabel('time (\mus)'); ylabel('amplitude'); title('Echo in white noise (one trial)');
subplot(2,1,2);
y1 = conv(X(:,1), flipud(s), 'same');
y0 = conv(x0, flipud(s), 'same');                % noise-free output
plot(tt, y1/sqrt(sigma^2*sum(s.^2)), 'k'); hold on;
plot(tt, y0/sqrt(sigma^2*sum(s.^2)), 'b', 'LineWidth', 1.2);
plot([0 tt(end)], sqrt(SNRth)*[1 1], 'r--');
legend('signal + noise', 'signal only', 'sqrt(2E/N0)', 'Location', 'southeast');
xlim([2 11]); xlabel('time (\mus)'); ylabel('output / noise std');
title(sprintf('Matched filter output: peak SNR theory %.1f dB', 10*log10(SNRth)));

figure('Position',[100 100 700 420]);
semilogx(r_list, 10*log10(effTh), 'b-', r_list, 10*log10(effEx), 'ko', r_list, 10*log10(effMC), 'r+', 'LineWidth', 1.2);
grid on; xlabel('filter bandwidth / echo bandwidth  r'); ylabel('SNR loss (dB)');
legend('2r/(1+r^2) narrowband theory', 'exact', 'Monte Carlo', 'Location', 'south');
title('Matched filters are flat near the optimum');

What to look for. \(\sum s^2/\sigma^2\) = 10.00 dB; Monte Carlo RF 9.94 dB, analytic \(\lvert y\rvert^2\) 6.95 dB (theory 6.99 = \(E/N_0\)), analytic \(\mathrm{Re}\{y\}\) 9.94 dB. Bandwidth ratio \(r\) = 0.5 and 2: −0.97 dB theory, −0.97 and −1.04 dB exact; \(r\) = 3: −2.22 theory vs −2.55 exact (broadband breakdown). Centre offsets 0.5, 1, 2 σ: −0.54, −2.17, −8.69 dB, equal to \(\exp(-D^2/2)\).

m10_chirp_compression.mLFM width and PSL vs TB against the closed form, taper table (receive, transmit, split), Monte Carlo compression gain · change B, TB_list, TBw, windows · runtime ≈ 3 s

%% m10_chirp_compression.m
% Linear FM (chirp) pulse compression: mainlobe width, range sidelobes,
% compression gain and the price of amplitude tapering.
%   - Rect-envelope LFM in complex baseband: the autocorrelation magnitude is
%     exactly (1-|tau|/T)|sinc(B tau (1-|tau|/T))|; the -6 dB width tends to
%     1.207/B and the peak sidelobe to -13.26 dB as TB grows.
%   - Compression gain at equal peak amplitude = energy ratio. Against a
%     Gaussian pulse with the same -6 dB bandwidth it is 1.505*TB.
%   - Tapers applied on receive only (mismatched filter), on transmit only,
%     or split as sqrt(w) on both sides.
% What to change: B, TB_list, TBw (taper table), the window list, Ntrial.
% Expected output: width*B 1.234 (TB 10) ... 1.208 (TB 200); PSL -14.9 ... -13.3 dB;
% Hann on receive at TB 50: PSL -31.3 dB, width*B 2.00, loss 1.76 dB; Hamming on
% receive: -36.5 dB (not -42.7: Fresnel ripple); Monte Carlo gain (10^4 trials) within about
% 0.25 dB of 10*log10(1.505*TB). Runtime about 3 s in Octave.
clear; close all; rng(1);

%% Helpers (anonymous only, so the script runs unchanged in MATLAB and Octave)
fs    = 200e6;                                    % complex-baseband sampling rate (Hz)
B     = 5e6;                                      % swept bandwidth (Hz)
tvec  = @(T) ((0:round(T*fs)-1)' - (round(T*fs)-1)/2)/fs;    % centred time axis (s)
lfm   = @(T) exp(1j*pi*(B/T)*tvec(T).^2);         % baseband up-chirp, instantaneous f = (B/T) t
sincf = @(x) (sin(pi*x) + (x==0))./(pi*x + (x==0));
nf    = @(n) 2^nextpow2(n);
% cross-correlation of received s with filter template h (output = s correlated with h), lag 0 centred
xcor  = @(s,h) fftshift(ifft(fft(s, nf(2*max(numel(s),numel(h)))).*conj(fft(h, nf(2*max(numel(s),numel(h)))))));
% -6 dB width (samples) with linear interpolation, around the peak index k0 of a normalized envelope a
xR    = @(a,i) i-1 + (a(i-1)-0.5)./(a(i-1)-a(i));
xL    = @(a,i) i+1 - (a(i+1)-0.5)./(a(i+1)-a(i));
w6f   = @(a,k0) xR(a, k0-1+find(a(k0:end) < 0.5, 1)) - xL(a, k0+1-find(a(k0:-1:1) < 0.5, 1));
% peak sidelobe (dB) outside the first minima on each side of the peak
pslf  = @(a,k0) 20*log10(max(a([1:(k0+1-find(diff(a(k0:-1:1)) > 0, 1))-1, (k0-1+find(diff(a(k0:end)) > 0, 1))+1:numel(a)])));
x01   = @(n) ((0:n-1)' + 0.5)/n;                  % window abscissa in (0,1)
wins  = {'rect','Tukey 25%','Hann','Hamming','Blackman'};
winf  = {@(x) ones(size(x)), ...
         @(x) 0.5 - 0.5*cos(2*pi*min(min(x,1-x),0.125)/0.25), ...   % Tukey: cosine ramps over 12.5% each side
         @(x) 0.5 - 0.5*cos(2*pi*x), ...
         @(x) 0.54 - 0.46*cos(2*pi*x), ...
         @(x) 0.42 - 0.5*cos(2*pi*x) + 0.08*cos(4*pi*x)};

%% 1. Rect LFM: width and PSL vs TB, against the closed form
TB_list = [10 20 50 100 200];
fprintf('Rect LFM, matched filter (complex baseband, B = %.0f MHz)\n   TB    width*B   PSL (dB)   closed-form max error\n', B/1e6);
for TBv = TB_list
    T = TBv/B; s = lfm(T);
    y = xcor(s,s); a = abs(y)/max(abs(y)); [~,k0] = max(a);
    tau = ((1:numel(a))' - k0)/fs;
    f = (abs(tau) < T).*(1 - abs(tau)/T).*abs(sincf(B*tau.*(1 - abs(tau)/T)));
    fprintf('%5d    %6.3f    %6.2f     %.1e\n', TBv, w6f(a,k0)/fs*B, pslf(a,k0), max(abs(a-f)));
end
fprintf('Limit TB -> inf: sinc, width*B = 1.207, PSL = -13.26 dB\n');

%% 2. Tapering: receive (mismatched), transmit, split; TB = 50
TBw = 50; T = TBw/B; s = lfm(T); n = numel(s); x = x01(n);
fprintf('\nTapers at TB = %d: PSL (dB) / width*B / SNR loss vs rect matched (dB)\n', TBw);
fprintf('  window       receive only            transmit only           split sqrt(w)\n');
ys = cell(1,numel(wins));
for k = 1:numel(wins)
    w = winf{k}(x); line = sprintf('  %-10s', wins{k});
    for mode = 1:3
        if mode == 1, st = s;           h = s.*w;        end     % mismatched receive filter
        if mode == 2, st = s.*w;        h = st;          end     % tapered transmit, matched receive
        if mode == 3, st = s.*sqrt(w);  h = st;          end     % split
        y = xcor(st,h); a = abs(y)/max(abs(y)); [~,k0] = max(a);
        loss = 10*log10(max(abs(y))^2/(sum(abs(h).^2)*n));     % rect matched reference: peak^2/sum|h|^2 = n
        line = [line sprintf('  %6.1f / %4.2f / %5.2f  ', pslf(a,k0), w6f(a,k0)/fs*B, loss)]; %#ok<AGROW>
        if mode == 1, ys{k} = a; end
    end
    fprintf('%s\n', line);
end
fprintf('Ideal spectral-window PSL for comparison: Hann -31.5, Hamming -42.7, Blackman -58.1 dB\n');

%% 3. Compression gain: Monte Carlo against a Gaussian pulse of the same bandwidth and peak amplitude
% A lower complex-baseband rate (5x the swept band) is enough here and keeps the Monte Carlo fast.
fs3 = 25e6;
lfm3 = @(T) exp(1j*pi*(B/T)*(((0:round(T*fs3)-1)' - (round(T*fs3)-1)/2)/fs3).^2);
sf = B/2.3548; st = 1/(2*pi*sf);                  % Gaussian amplitude spectrum std (Hz), envelope std (s)
tg = (-ceil(5*st*fs3):ceil(5*st*fs3))'/fs3; g = exp(-tg.^2/(2*st^2));   % reference pulse, peak 1
Ntrial = 10000; sigma = 3;                        % complex white noise, E|n|^2 = sigma^2
fprintf('\nCompression gain (equal peak amplitude 1, same -6 dB bandwidth)\n   TB   theory 1.505*TB   energy ratio   Monte Carlo\n');
for TBv = [10 50 100]
    T = TBv/B; s = lfm3(T);
    SNR = @(v) abs(mean(v))^2/mean(abs(v-mean(v)).^2);
    nz = sigma*(randn(numel(s),Ntrial) + 1j*randn(numel(s),Ntrial))/sqrt(2);
    yc = s'*(repmat(s,1,Ntrial) + nz);            % matched filter output at the peak
    nz = sigma*(randn(numel(g),Ntrial) + 1j*randn(numel(g),Ntrial))/sqrt(2);
    yg = g'*(repmat(g,1,Ntrial) + nz);
    fprintf('%5d     %6.2f dB        %6.2f dB      %6.2f dB\n', TBv, 10*log10(1.505*TBv), ...
        10*log10(sum(abs(s).^2)/sum(g.^2)), 10*log10(SNR(yc)/SNR(yg)));
end

%% Figures
T = 10e-6; s = lfm(T); n = numel(s); tt = tvec(T)*1e6;
figure('Position',[100 100 900 600]);
subplot(2,1,1);
f0 = 5e6; tr = tvec(T); plot(tr*1e6, real(s.*exp(1j*2*pi*f0*tr)), 'b'); xlim([-5 5]);
xlabel('time (\mus)'); ylabel('amplitude');
title(sprintf('RF chirp, f_0 = 5 MHz, B = %.0f MHz, T = %.0f \\mus (TB = %d)', B/1e6, T*1e6, round(T*B)));
subplot(2,1,2);
N = 2^16; S = fftshift(abs(fft(s, N))); fax = ((0:N-1)' - N/2)*fs/N;
plot(fax/1e6 + 5, S/max(S), 'b'); xlim([1 9]); grid on;
xlabel('frequency (MHz)'); ylabel('|S(f)| (normalized)');
title('Spectrum: nearly rectangular, Fresnel ripple at the band edges');

figure('Position',[100 100 900 420]);
lag = ((1:numel(ys{1}))' - (numel(ys{1})/2+1))/fs*1e6;
plot(lag, 20*log10(ys{1}+1e-12), 'k', lag, 20*log10(ys{3}+1e-12), 'b', lag, 20*log10(ys{4}+1e-12), 'r');
xlim([-10 10]); ylim([-80 3]); grid on;
xlabel('lag (\mus)'); ylabel('envelope (dB)');
legend('rect (matched)', 'Hann on receive', 'Hamming on receive', 'Location', 'northeast');
title(sprintf('Compressed chirp, TB = %d', TBw));

What to look for. width·B 1.234, 1.221, 1.212, 1.210, 1.208 and PSL −14.90, −13.96, −13.52, −13.39, −13.33 dB for TB 10–200, closed-form error ≤ \(3\times10^{-5}\). The taper table of 10.3. Compression gain at TB 10, 50, 100: theory 11.78, 18.77, 21.78 dB; energy ratio identical; Monte Carlo 11.86, 18.70, 21.56 dB.

m10_barker_golay.mBarker-13 matched and LS mismatched, Golay pairs N = 4–64, motion between firings (phase model and RF simulation) · change Lmm, Nlist, f0, fbw, d_um · runtime ≈ 1 s

%% m10_barker_golay.m
% Binary phase codes: Barker-13 (matched and least-squares mismatched filter)
% and Golay complementary pairs, including what motion between the two
% firings does to the Golay sidelobe cancellation.
%   - Barker-13: autocorrelation peak 13, sidelobes 0 or 1 -> PSL = -22.28 dB.
%   - LS mismatched filter of length L: minimizes the summed squared sidelobes;
%     PSL falls with L for about 0.2 dB of SNR loss.
%   - Golay pair (A,B): R_A + R_B = 2N delta. If echo B is phase-rotated by phi
%     (axial displacement d between firings, phi = 4 pi f0 d / c), the residual
%     PSL is (max|R_A(k)|/N) tan(phi/2).
%   - RF check: 1-cycle chips at f0, Gaussian transducer, true delay 2d/c.
% What to change: Lmm, Nlist, f0, fbw, d_um.
% Expected output: Barker PSL -22.28 dB, gain 11.14 dB; LS filter L = 39: PSL
% -38.5 dB, loss 0.20 dB; Golay N = 16: single-code PSL -10.1 dB, pair sum
% sidelobes exactly 0, d = 2 um -> -37.9 dB (phase-only theory), -38.4 dB on RF.
% Runtime about 1 s.
clear; close all; rng(1);

dB = @(x) 20*log10(max(abs(x), 1e-300));
c  = 1540;                 % sound speed (m/s)
f0 = 5e6;                  % centre frequency (Hz)

%% 1. Barker-13, matched filter
b13 = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1]';
r = conv(b13, flipud(b13));                      % autocorrelation, peak at index 13
side = r; side(13) = [];
fprintf('Barker-13 autocorrelation: %s\n', mat2str(r'));
fprintf('PSL = %.2f dB, compression gain 10log10(13) = %.2f dB\n', dB(max(abs(side))/13), 10*log10(13));

%% 2. Least-squares mismatched filters for Barker-13
Lmm = [13 21 27 39 51 65];
fprintf('\nLS mismatched filter (target: unit peak, zero sidelobes)\n    L    PSL (dB)   SNR loss (dB)\n');
hbest = [];
for L = Lmm
    M = 13 + L - 1; C = zeros(M, L);
    for j = 1:L, C(j:j+12, j) = b13; end           % convolution matrix: y = C*h
    d = zeros(M,1); d((M+1)/2) = 1;                % desired output: a centred spike
    h = C \ d; y = C*h;
    [pk, k] = max(abs(y)); ys = y; ys(k) = [];
    loss = 10*log10(pk^2/(sum(h.^2)*sum(b13.^2)));  % relative to the matched filter (13)
    fprintf('%5d   %7.1f     %6.2f\n', L, dB(max(abs(ys))/pk), loss);
    if L == 39, hbest = h; ybest = y/pk; end
end

%% 3. Golay complementary pairs
Nlist = [4 8 16 32 64];
fprintf('\nGolay pairs\n    N   max|R_A|   PSL single (dB)   max|R_A+R_B| off-peak   gain 10log10(2N)\n');
for N = Nlist
    a = 1; b = 1;
    while numel(a) < N, t = [a b]; b = [a -b]; a = t; end   % recursive construction
    ra = conv(a, fliplr(a)); rb = conv(b, fliplr(b)); s = ra + rb;
    sa = ra; sa(N) = []; ss = s; ss(N) = [];
    fprintf('%5d   %5d       %7.1f             %g                 %.2f dB\n', N, max(abs(sa)), dB(max(abs(sa))/N), max(abs(ss)), 10*log10(2*N));
end
N = 16; a = 1; b = 1;
while numel(a) < N, t = [a b]; b = [a -b]; a = t; end
ra = conv(a, fliplr(a)); rb = conv(b, fliplr(b)); sa = ra; sa(N) = []; mA = max(abs(sa));
d_um = [0.5 1 2 5 10 20 38.5];
fprintf('\nGolay N = 16, phase rotation phi = 4 pi f0 d / c between firings (f0 = %.0f MHz)\n   d (um)   phi (rad)   PSL simulated   PSL theory (mA/N) tan(phi/2)\n', f0/1e6);
for d = d_um
    phi = 4*pi*f0*d*1e-6/c; z = ra + exp(1j*phi)*rb; zs = z; zs(N) = [];
    fprintf('%8.1f   %7.3f      %7.1f dB        %7.1f dB\n', d, phi, dB(max(abs(zs))/abs(z(N))), dB(mA/N*tan(phi/2)));
end

%% 4. RF simulation: 1-cycle chips, transducer band limit, true delay 2d/c on firing B
fs = 16*f0; nc = round(fs/f0);                    % 16 samples per chip (one carrier cycle)
fbw = 0.7; sf = fbw*f0/2.3548;                    % two-way Gaussian transducer response
chip = sin(2*pi*(0:nc-1)'/nc);
codewave = @(code) reshape(chip*code(:)', [], 1);
Nfft = 4096; fax = [0:Nfft/2, -Nfft/2+1:-1]'*fs/Nfft;
Htr = exp(-(abs(fax)-f0).^2/(2*sf^2));
wA = real(ifft(fft(codewave(a), Nfft).*Htr)); wB = real(ifft(fft(codewave(b), Nfft).*Htr));   % received pulse shapes
an = zeros(Nfft,1); an(1) = 1; an(2:Nfft/2) = 2; an(Nfft/2+1) = 1;     % analytic-signal weights
env = @(Y) abs(ifft(Y.*an));
WA = fft(wA); WB = fft(wB);
fprintf('\nRF simulation (N = 16, 1 cycle/chip, %.0f%% transducer)\n   d (um)   PSL sum (dB)   phase-only theory (dB)\n', fbw*100);
dplot = [0 5 20]; envs = cell(1,3);
for d = [0 0.5 2 5 20]
    tau = 2*d*1e-6/c;                             % extra round-trip delay of echo B (s)
    YB = WB.*exp(-1j*2*pi*fax*tau);               % exact fractional delay (band-limited signal)
    Ysum = WA.*conj(WA) + YB.*conj(WB);           % each firing compressed by its own matched filter
    e = fftshift(env(Ysum)); [pk,k] = max(e); en = e/pk;
    % sidelobe region: more than 4 chips from the peak (the band-limited mainlobe skirt is excluded)
    mask = abs((1:Nfft)' - k) > 4*nc;
    psl = dB(max(en(mask)));
    phi = 4*pi*f0*d*1e-6/c;
    if d == 0, th = -Inf; else, th = dB(mA/N*tan(phi/2)); end
    fprintf('%8.1f      %7.1f          %7.1f\n', d, psl, th);
    ii = find(dplot == d); if ~isempty(ii), envs{ii} = en; end
end
eA = fftshift(env(WA.*conj(WA))); eA = eA/max(eA);

%% Figures
figure('Position',[100 100 900 560]);
subplot(2,1,1);
lag = (-12:12)';
plot(lag, max(dB(r/13), -80), 'ko-', 'LineWidth', 1.2); hold on;   % zeros drawn at -80 dB
lagm = (1:numel(ybest))' - (numel(ybest)+1)/2;
plot(lagm, max(dB(ybest), -80), 'r.-');
ylim([-80 3]); grid on; xlabel('lag (chips)'); ylabel('output (dB)');
legend('Barker-13 matched (-22.3 dB)', 'LS mismatched, L = 39', 'Location', 'northeast');
title('Barker-13: matched vs least-squares mismatched filter');
subplot(2,1,2);
tt = ((1:Nfft)' - Nfft/2 - 1)/fs*1e6;
plot(tt, dB(eA), 'Color', [0.6 0.6 0.6]); hold on;
plot(tt, dB(envs{1}), 'k', tt, dB(envs{2}), 'b', tt, dB(envs{3}), 'r');
xlim([-4 4]); ylim([-70 3]); grid on; xlabel('lag (\mus)'); ylabel('envelope (dB)');
legend('A alone', 'A+B static', 'A+B, d = 5 \mum', 'A+B, d = 20 \mum', 'Location', 'northeast');
title('Golay N = 16 on RF: motion between firings breaks the cancellation');

What to look for. Barker PSL −22.28 dB, gain 11.14 dB; LS filter \(L\) = 13, 39, 65: −24.0, −38.5, −58.6 dB for 0.14–0.21 dB. Golay \(N\) = 4…64: single-code PSL −12.0, −8.5, −10.1, −13.2, −13.8 dB and pair-sum sidelobes exactly 0. \(N\) = 16 motion: 2 µm −37.9, 5 µm −29.9, 20 µm −17.4 dB (simulation = theory); RF with 70 % transducer: −38.4, −30.4, −17.9 dB, static −102 dB.

m10_attenuation_mismatch.mChirp × transducer × attenuation, fixed vs depth-adapted vs equalizing filters, centroid shift of the raw echo · change alpha_dB, zlist, fbw, T, B, epsW · runtime ≈ 1 s

%% m10_attenuation_mismatch.m
% A chirp through attenuating tissue, compressed with (a) the filter matched to
% the transmitted echo, (b) a depth-adapted matched filter, (c) an equalizing
% (Wiener-type) filter that restores the surface pulse shape.
%   - Attenuation alpha (dB/cm/MHz), linear in f, two-way path 2z, no dispersion.
%   - Received spectrum E(f) = S(f) H_tr(f) exp(-2 alpha_Np f z).
%   - SNR loss of (a) vs (b) = (int|E0||E|)^2 / (int|E0|^2 int|E|^2) <= 1 (Cauchy-Schwarz).
%   - With an exact (conjugate) filter the output spectrum is real and >= 0,
%     so the compressed envelope does not move; the raw echo's energy centroid
%     moves by about (fc - f0) T / B.
% What to change: alpha_dB, zlist, fbw, T, B, epsW.
% Expected output (alpha 0.5, 80 % transducer): centroid 5.00 -> 3.68 MHz at 6 cm,
% mismatch loss 2.0 dB at 6 cm and 4.6 dB at 10 cm; envelope shift 0.0 ns at
% every depth; at 6 cm the equalizer holds the -6 dB width at 0.38 us (0.34 at the
% surface, 0.48 for the adapted filter) but costs 5.8 dB.
% The adapted filter's PSL looks good only because its broad skirt has no null:
% read the -20 dB width (w20) as well. Runtime about 1 s.
clear; close all; rng(1);

%% Parameters
fs  = 50e6;                    % sampling rate (Hz)
f0  = 5e6;                     % chirp centre frequency (Hz)
T   = 10e-6;                   % chirp duration (s)
B   = 5e6;                     % swept bandwidth (Hz), 2.5-7.5 MHz
fbw = 0.8;                     % -6 dB fractional bandwidth of the two-way transducer response
alpha_dB = 0.5;                % attenuation (dB/cm/MHz)
zlist = [0 2 4 6 8 10];        % depth (cm)
epsW  = 1e-2;                  % equalizer regularization, relative to max|E|^2 (-20 dB)

Nfft = 4096; fax = [0:Nfft/2, -Nfft/2+1:-1]'*fs/Nfft;
n  = round(T*fs); t = ((0:n-1)' - (n-1)/2)/fs;
s  = cos(2*pi*(f0*t + B/(2*T)*t.^2));                 % RF up-chirp
S  = fft(s, Nfft);
sf = fbw*f0/2.3548;
Htr = exp(-(abs(fax)-f0).^2/(2*sf^2));                % Gaussian two-way transducer response
E0 = S.*Htr;                                          % echo spectrum at z = 0
an = zeros(Nfft,1); an(1) = 1; an(2:Nfft/2) = 2; an(Nfft/2+1) = 1;
env = @(Y) fftshift(abs(ifft(Y.*an)));                % envelope of a filter output, lag 0 centred
k0 = Nfft/2 + 1;
xR  = @(a,i,v) i-1 + (a(i-1)-v)./(a(i-1)-a(i));
xL  = @(a,i,v) i+1 - (a(i+1)-v)./(a(i+1)-a(i));
wf  = @(a,k,v) xR(a, k-1+find(a(k:end) < v, 1), v) - xL(a, k+1-find(a(k:-1:1) < v, 1), v);  % full width at level v
pslf = @(a,k) 20*log10(max(a([1:(k+1-find(diff(a(k:-1:1)) > 0, 1))-1, (k-1+find(diff(a(k:end)) > 0, 1))+1:numel(a)])));
pos = fax > 0; tt = (0:Nfft-1)'/fs;

fprintf('alpha = %.2f dB/cm/MHz, chirp %.1f-%.1f MHz, T = %.0f us, transducer %.0f%%\n', alpha_dB, (f0-B/2)/1e6, (f0+B/2)/1e6, T*1e6, fbw*100);
fprintf(' z(cm)  fc(MHz)  2*a*f0*z(dB) | fixed: w6 w20 (ns) PSL(dB) loss(dB) shift(ns) | adapted: w6 w20 PSL | equalizer: w6 w20 PSL loss | echo centroid shift (us) vs (fc-f0)T/B\n');
keep = cell(3,1);
for z = zlist
    E = E0.*exp(-alpha_dB/8.686*abs(fax)/1e6*2*z); % amplitude attenuation, alpha in Np/cm/MHz, path 2z
    P = abs(E).^2; fc = sum(fax(pos).*P(pos))/sum(P(pos));
    Hs = {E0, E, E.*abs(E0).^2./(abs(E).^2 + epsW*max(abs(E).^2))};
    w6 = zeros(1,3); w20 = w6; psl = w6; snr = w6; sh = w6;
    for k = 1:3
        e = env(E.*conj(Hs{k})); [pk, kp] = max(e); a = e/pk;
        w6(k) = wf(a,kp,0.5)/fs; w20(k) = wf(a,kp,0.1)/fs; psl(k) = pslf(a,kp); sh(k) = (kp - k0)/fs;
        snr(k) = pk^2/sum(abs(Hs{k}).^2);           % peak^2 / output noise power (common scale)
        if z == 6, keep{k} = a; end
    end
    ea = abs(ifft(E.*an)).^2; ea0 = abs(ifft(E0.*an)).^2;    % raw echo envelopes (no compression)
    dtc = sum(tt.*ea)/sum(ea) - sum(tt.*ea0)/sum(ea0);
    fprintf(' %4.0f   %5.2f     %5.1f      |  %4.0f %4.0f  %6.1f  %6.2f   %5.2f   | %4.0f %5.0f %6.1f | %4.0f %4.0f %6.1f %6.2f |  %6.2f vs %6.2f\n', ...
        z, fc/1e6, 2*alpha_dB*f0/1e6*z, w6(1)*1e9, w20(1)*1e9, psl(1), 10*log10(snr(1)/snr(2)), sh(1)*1e9, ...
        w6(2)*1e9, w20(2)*1e9, psl(2), w6(3)*1e9, w20(3)*1e9, psl(3), 10*log10(snr(3)/snr(2)), dtc*1e6, (fc-f0)*T/B*1e6);
    if z == 6, E6 = E; end
end
alpha0 = alpha_dB/8.686;       % Np/cm/MHz
fprintf('Gaussian-spectrum shift for reference: fc = f0 - 2 alpha0 sigma_f^2 z, sigma_f = %.2f MHz -> %.2f MHz at 6 cm (the flat chirp band makes the true shift smaller)\n', ...
    sf/1e6, f0/1e6 - 2*alpha0*(sf/1e6)^2*6);

%% Figures
figure('Position',[100 100 900 600]);
subplot(2,1,1);
plot(fax(pos)/1e6, abs(E0(pos))/max(abs(E0)), 'k', fax(pos)/1e6, abs(E6(pos))/max(abs(E6)), 'r', 'LineWidth', 1.2);
xlim([0 10]); grid on; xlabel('frequency (MHz)'); ylabel('|E(f)| (normalized)');
legend('z = 0', 'z = 6 cm', 'Location', 'northeast');
title(sprintf('Echo spectrum: chirp x transducer x attenuation (%.1f dB/cm/MHz)', alpha_dB));
subplot(2,1,2);
lag = ((1:Nfft)' - k0)/fs*1e6;
plot(lag, 20*log10(keep{1}+1e-9), 'k', lag, 20*log10(keep{2}+1e-9), 'b', lag, 20*log10(keep{3}+1e-9), 'r');
xlim([-12 12]); ylim([-70 3]); grid on; xlabel('lag (\mus)'); ylabel('envelope (dB)');
legend('fixed (matched to z = 0)', 'depth-adapted matched', 'equalizer to z = 0 shape', 'Location', 'northeast');
title('Compressed echo at 6 cm');

What to look for. At 0, 2, 4, 6, 8, 10 cm: centroid 5.00, 4.49, 4.04, 3.68, 3.40, 3.16 MHz; fixed-filter loss 0, −0.26, −0.97, −2.01, −3.23, −4.55 dB and envelope shift 0.00 ns throughout. At 6 cm: −6 dB widths 374 (fixed), 475 (adapted), 376 ns (equalizer, −5.77 dB); PSL −26.1, −35.8 (no null), −29.0 dB. Raw-echo centroid shift −2.47 µs vs \((f_c - f_0)T/B\) = −2.64 µs (the code length truncates it at depth).

m10_snr_budget_output.mPenetration gain table, broadband depth check, MI / ISPTA.3 / TI at equal peak pressure, equal-energy averaging · change alpha_dB, f0, fbw, TB, PRF · runtime ≈ 1 s

%% m10_snr_budget_output.m
% SNR budget of coded excitation and its acoustic-output cost.
%   - Penetration gain: Delta z = G / (2 alpha f0) with alpha in dB/cm/MHz,
%     two-way path, narrowband. Broadband check: SNR(z) = int |E0(f)|^2 exp(-4 alpha0 f z) df.
%   - Acoustic output at equal peak pressure: MI is unchanged; PII, ISPTA and
%     time-averaged power (hence TI) rise by the energy ratio, which is the
%     coding gain itself.
%   - At equal total energy, one coded firing and K averaged short pulses give
%     the same SNR (white noise): coding saves time and peak pressure, not energy.
% What to change: alpha_dB, f0, fbw, TB, PRF, p_r3 (derated peak pressure).
% Expected output: 15 dB at 0.5 dB/cm/MHz and 5 MHz -> 3.0 cm; a broadband
% (Gaussian, cut at -24 dB) spectrum gives 4.9 cm because it downshifts (2.8 MHz there);
% chirp TB 50 at MI 1.9 raises ISPTA.3 at a point (200 pulses/s) from 27 to
% 2003 mW/cm^2 (x75.2 = 18.8 dB, the coding gain; limit 720); one coded firing
% vs 75 averaged short pulses: SNR equal within Monte Carlo error (about 0.2 dB).
% Runtime about 3 s.
clear; close all; rng(1);

c = 1540;            % sound speed (m/s)
rhoc = 1.5e6;        % characteristic impedance of tissue (Rayl), for intensity p^2/(rho c)

%% 1. Penetration gain, narrowband rule
codes = {'chirp TB 50, rect', 'chirp TB 50, Hann on receive', 'Barker-13', 'Golay 16 (pair, 2 firings)'};
G = [10*log10(1.505*50), 10*log10(1.505*50) - 1.76, 10*log10(13), 10*log10(2*16)];
fprintf('Penetration gain Delta z = G/(2 alpha f0)  (cm)\n%-30s  G (dB)   a=0.5,f0=3   a=0.5,f0=5   a=0.7,f0=5   a=0.5,f0=7.5\n', 'code');
for k = 1:numel(codes)
    fprintf('%-30s  %5.2f     %6.2f       %6.2f       %6.2f       %6.2f\n', codes{k}, G(k), ...
        G(k)/(2*0.5*3), G(k)/(2*0.5*5), G(k)/(2*0.7*5), G(k)/(2*0.5*7.5));
end
fprintf('Reference: 15 dB at 0.5 dB/cm/MHz, 5 MHz -> %.2f cm\n', 15/(2*0.5*5));

%% 2. Broadband check: the attenuation slope flattens with depth because the spectrum downshifts
alpha_dB = 0.5; f0 = 5; fbw = 0.6;                    % dB/cm/MHz, MHz, -6 dB fractional bandwidth
sf = fbw*f0/2.3548; f = linspace(0.01, 12, 1200)';   % MHz
E0sq = exp(-(f-f0).^2/sf^2).*(abs(f-f0) <= fbw*f0);  % |E0(f)|^2: Gaussian amplitude spectrum cut at its -24 dB points
a0 = alpha_dB/8.686;                                  % Np/cm/MHz
z = linspace(0, 30, 601);
SNRdB = 10*log10(sum(repmat(E0sq,1,numel(z)).*exp(-4*a0*f*z), 1)/sum(E0sq));   % relative to z = 0
SNRnb = -2*alpha_dB*f0*z;                             % narrowband rule
budget = 40;                                          % short pulse: detectable down to -40 dB re z = 0
zsp_nb = budget/(2*alpha_dB*f0); zsp_bb = interp1(-SNRdB, z, budget);
fprintf('\nBroadband check (alpha %.1f, f0 %.0f MHz, %.0f%% band), short-pulse budget %d dB\n', alpha_dB, f0, fbw*100, budget);
fprintf('   G (dB)   depth gain narrowband   depth gain broadband   centre frequency there (MHz)\n');
for Gk = [10 15 20]
    zG = interp1(-SNRdB, z, budget+Gk);
    fprintf('   %4.0f        %5.2f cm              %5.2f cm               %4.2f\n', Gk, Gk/(2*alpha_dB*f0), zG - zsp_bb, sum(f.*E0sq.*exp(-4*a0*f*zG))/sum(E0sq.*exp(-4*a0*f*zG)));
end
fprintf('   short-pulse depth: narrowband %.2f cm, broadband %.2f cm\n', zsp_nb, zsp_bb);

%% 3. Acoustic output at equal derated peak pressure
fc  = 5e6; B = fbw*fc;                                % Hz
p_r3 = 1.9*sqrt(fc/1e6)*1e6;                          % derated peak rarefactional pressure at MI = 1.9 (Pa)
fs  = 200e6;
PRF = 200;           % effective pulse rate at one point (Hz): scanned B-mode, e.g. 40 frames/s x 5 overlapping lines
sft = B/2.3548; st = 1/(2*pi*sft);
tg  = (-6*st:1/fs:6*st)'; pg = -p_r3*exp(-tg.^2/(2*st^2)).*cos(2*pi*fc*tg);     % short pulse, rarefactional peak p_r3
T   = 50/B; tc = (-T/2:1/fs:T/2)'; pc = p_r3*cos(2*pi*(fc*tc + B/(2*T)*tc.^2)); % chirp TB 50
PII = @(p) sum(p.^2)/fs/rhoc;                         % pulse-intensity integral (J/m^2)
MI  = @(p) max(-p)/1e6/sqrt(fc/1e6);
Ig = PRF*PII(pg)/10; Ic = PRF*PII(pc)/10;             % W/m^2 -> mW/cm^2 is x0.1
fprintf('\nAcoustic output at one point, effective pulse rate %.0f Hz, derated peak %.2f MPa (MI 1.9 at %.0f MHz)\n', PRF, p_r3/1e6, fc/1e6);
fprintf('   short pulse: MI %.2f, ISPTA.3 %.0f mW/cm^2\n', MI(pg), Ig);
fprintf('   chirp TB 50: MI %.2f, ISPTA.3 %.0f mW/cm^2  (ratio %.1f = %.2f dB; coding gain %.2f dB)\n', MI(pc), Ic, Ic/Ig, 10*log10(Ic/Ig), 10*log10(1.505*50));
fprintf('   to meet ISPTA.3 <= 720 mW/cm^2: pulse rate <= %.0f Hz, or peak pressure x %.2f (MI %.2f)\n', PRF*720/Ic, min(1,sqrt(720/Ic)), MI(pc)*min(1,sqrt(720/Ic)));
fprintf('   same pulses at 5 kHz on every point (unfocused ultrafast-like, same peak): short %.0f, chirp %.0f mW/cm^2\n', Ig*5e3/PRF, Ic*5e3/PRF);
fprintf('   TI scales with time-averaged power: x %.1f at the same PRF\n', Ic/Ig);
fprintf('   max PRF set by depth: 5 cm %.0f Hz, 10 cm %.0f Hz, 15 cm %.0f Hz\n', c/(2*0.05), c/(2*0.10), c/(2*0.15));
for Tk = [5 10 20 40]
    fprintf('   blind zone c T / 2 for T = %2d us: %.1f mm\n', Tk, c*Tk*1e-6/2*1e3);
end

%% 4. Equal energy, equal SNR: one coded firing vs K averaged short pulses
K = round(1.505*50);                                  % energy ratio of the TB-50 chirp to the short pulse
Ntrial = 2000; sigma = 4;                             % complex white noise
fs = 25e6;                                            % complex baseband rate is enough here (B = 3 MHz)
tb = (-6*st:1/fs:6*st)'; g = exp(-tb.^2/(2*st^2));    % baseband short pulse, peak 1
tcb = (-T/2:1/fs:T/2)'; sc = exp(1j*pi*(B/T)*tcb.^2); % baseband chirp, peak 1
SNR = @(v) abs(mean(v))^2/mean(abs(v-mean(v)).^2);
yc = sc'*(repmat(sc,1,Ntrial) + sigma*(randn(numel(sc),Ntrial)+1j*randn(numel(sc),Ntrial))/sqrt(2));
yg = zeros(1,Ntrial);
for k = 1:K
    yg = yg + g'*(repmat(g,1,Ntrial) + sigma*(randn(numel(g),Ntrial)+1j*randn(numel(g),Ntrial))/sqrt(2));
end
fprintf('\nEqual energy: coded single firing %.2f dB, %d averaged short pulses %.2f dB (energy ratio %.2f dB)\n', ...
    10*log10(SNR(yc)), K, 10*log10(SNR(yg)), 10*log10(sum(abs(sc).^2)/sum(g.^2)));

%% Figures
figure('Position',[100 100 900 420]);
plot(z, SNRnb, 'k--', z, SNRdB, 'b', z, SNRdB + 15, 'r', 'LineWidth', 1.2); hold on;
plot([0 20], -budget*[1 1], 'Color', [0.5 0.5 0.5]);
ylim([-120 20]); grid on; xlabel('depth (cm)'); ylabel('echo SNR re short pulse at z = 0 (dB)');
legend('narrowband -2\alpha f_0 z', 'broadband short pulse', 'coded, +15 dB', 'detection threshold', 'Location', 'southwest');
title(sprintf('SNR budget: %.1f dB/cm/MHz, %.0f MHz, %.0f%% band', alpha_dB, f0, fbw*100));

What to look for. 15 dB → 3.00 cm at 0.5 dB/cm/MHz and 5 MHz; TB-50 rect 3.75 cm (6.26 cm at 3 MHz). Broadband, 40 dB budget: short-pulse depth 9.66 cm (narrowband 8.00), gains of 10/15/20 dB add 3.16/4.88/6.70 cm with centre frequencies 2.98/2.83/2.70 MHz. Output: MI 1.90 for both, \(I_{\mathrm{SPTA.3}}\) 27 → 2003 mW/cm² (ratio 75.2 = 18.76 dB), pulse rate ≤ 72 Hz or peak × 0.60 (MI 1.14) to meet 720 mW/cm². Equal energy: 14.36 vs 14.12 dB.

10.13Self-check

Q1A real RF echo has Σs² = 100 (sampled) in white noise of variance σ² = 4. What is the matched-filter peak SNR, and what does a complex |y|² definition report on the analytic signal?

\(\sum s^2/\sigma^2\) = 25 = 14.0 dB (\(= 2E/N_0\)). The analytic \(\lvert y\rvert^2/\mathrm{E}\lvert n\rvert^2\) definition reports 12.5 = 11.0 dB (\(E/N_0\)). Same data, same detectability at high SNR.

Q2A chirp has T = 20 µs and B = 4 MHz. Give TB, the compressed −6 dB width and axial resolution, the gain against a Gaussian pulse of the same band, and the blind zone.

\(TB\) = 80. Width \(1.207/4\ \mathrm{MHz}\) = 302 ns, \(\Delta z = c \cdot 302\ \mathrm{ns}/2\) = 0.232 mm. \(G = 10\log_{10}(1.505 \times 80)\) = 20.8 dB. Blind zone \(cT/2\) = 15.4 mm.

Q3What does a Hann receive taper do to a TB-50 chirp?

PSL −13.5 → −31.3 dB, −6 dB width \(1.21/B \to 2.00/B\) (×1.65), SNR loss \((\sum w)^2/(N\sum w^2)\) = 0.667 = −1.76 dB.

Q4A Golay pair with N = 16 at 5 MHz sees 5 µm of axial motion between firings. What is the residual PSL, and what velocity is that at a 5 kHz PRF?

\(\varphi = 4\pi \cdot 5\ \mathrm{MHz} \cdot 5\ \mu\mathrm{m}/1540\ \mathrm{m/s}\) = 0.204 rad; \(\mathrm{PSL} = (5/16)\tan(0.102)\) = 0.032 → −29.9 dB. Velocity 5 µm × 5 kHz = 25 mm/s.

Q5How much depth does an 18.8 dB compression gain buy at 0.5 dB/cm/MHz at 3 MHz and at 7.5 MHz?

\(18.77/(2 \times 0.5 \times 3)\) = 6.3 cm; \(18.77/(2 \times 0.5 \times 7.5)\) = 2.5 cm (narrowband rule, noise-limited region).

Q6MI 1.9 at 5 MHz, 200 pulses per second at a point. Compare I_SPTA.3 for a 60 % Gaussian pulse and a TB-50 chirp of the same band and peak pressure.

\(p_{r.3} = 1.9\sqrt5\) = 4.25 MPa. Pulse: \(\mathrm{PII} = p^2\sqrt\pi\,\sigma_t/(2\rho c)\) = 1.33 J/m² → 27 mW/cm². Chirp: \(\mathrm{PII} = p^2T/(2\rho c)\) with \(T\) = 16.7 µs = 100 J/m² → about 2000 mW/cm², ratio 75 = \(1.505\,TB\), above the 720 mW/cm² limit. MI is the same for both.

Q7Why does the compressed chirp peak stay put under attenuation while the raw echo moves 2.5 µs earlier at 6 cm?

Linear-in-f attenuation without dispersion is a real, non-negative spectral weighting, so \(E_zE_0^*\) is real and ≥ 0 and the compressed envelope peaks at lag 0. The raw up-chirp emits low frequencies first; weighting them up moves its energy centroid earlier by \(\approx (f_c - f_0)T/B\).

Q8Why can a biphase code not be used directly for second-harmonic imaging, while a chirp can?

The second harmonic carries twice the phase: 0 and \(\pi\) both map to 0 (mod \(2\pi\)), so the code disappears and the harmonic is an uncompressible long pulse. A chirp's doubled phase is still quadratic: a chirp at \(2f_0\) sweeping \(2B\), compressible with its own matched filter.

10.14Figures from the literature

Published figures with measured data that show compression gain, sidelobes and the equal-energy argument outside the idealized models here. Links only; no figures are reproduced.

Measured axial PSF of short pulse, chirp and Barker. Fig. 8 shows measured axial point-spread functions and lateral speckle functions at the focal depth for a short pulse, a chirp and Barker codes on a real scanner; Fig. 5 compares short pulse, chirp, 7-bit Barker and 8-bit Golay transfer functions, and Fig. 9 the resulting strain images of a lesion phantom. Look at the range sidelobe floor of the measured chirp against the −13 dB rect-chirp prediction of (10.9): the transducer has already tapered it.

J. Liu, M. F. Insana, "Coded pulse excitation for ultrasonic strain imaging," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52(2):231–240, 2005, Figs. 5, 8 and 9. doi:10.1109/TUFFC.2005.1406549

Chirp versus short pulse on a phantom, SNR against depth. Fig. 13 shows tissue-mimicking phantom images with short-pulse and chirp-coded excitation; Fig. 14 plots the measured SNR as a function of depth for both; Figs. 5–7 show Tukey and Gaussian windowed chirps and their spectra before and after the matched filter. Compare the depth gap between the curves with \(\Delta z = G/(2\alpha f_0)\), (10.20). Open access (CC BY 4.0).

G.-C. Chun, H.-J. Chiang, K.-H. Lin, C.-M. Li, P.-J. Chen, T. Chen, "Ultrasound elasticity imaging system with chirp-coded excitation for assessing biomechanical properties of elasticity phantom," Materials 8(12):8392–8413, 2015, Figs. 5–7, 13 and 14. doi:10.3390/ma8125458

Equal energy, equal SNR, with spatial (Hadamard) coding. Figs. 3(a) and 4(a) give measured SNR versus depth in a calibrated phantom for multiplane-wave (Hadamard-coded) transmission against conventional plane-wave compounding with the same number of firings; Fig. 7 shows in vivo ultrafast Doppler of a rat brain. The gain comes from delivering more energy per firing at the same peak pressure, the argument of 10.7 applied across angles instead of time. Open access (CC BY 3.0).

E. Tiran et al., "Multiplane wave imaging increases signal-to-noise ratio in ultrafast ultrasound imaging," Phys. Med. Biol. 60:8549, 2015, Figs. 3, 4 and 7. doi:10.1088/0031-9155/60/21/8549

10.15Go deeper

G. L. Turin, “An introduction to matched filters,” IRE Trans. Inf. Theory, 6(3):311–329, 1960. The classic derivation, including colored noise.
J. R. Klauder, A. C. Price, S. Darlington, W. J. Albersheim, “The theory and design of chirp radars,” Bell Syst. Tech. J., 39(4):745–808, 1960. LFM, stationary phase, Fresnel ripple, weighting.
M. J. E. Golay, “Complementary series,” IRE Trans. Inf. Theory, 7(2):82–87, 1961.
R. H. Barker, “Group synchronizing of binary digital systems,” in W. Jackson (ed.), Communication Theory, Butterworths, 1953.
M. O'Donnell, “Coded excitation system for improving the penetration of real-time phased-array imaging systems,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 39(3):341–351, 1992.
B. Haider, P. A. Lewin, K. E. Thomenius, “Pulse elongation and deconvolution filtering for medical ultrasonic imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 45(1):98–113, 1998.
T. Misaridis, J. A. Jensen, “Use of modulated excitation signals in medical ultrasound. Part I: Basic concepts and expected benefits,” 52(2):177 (doi:10.1109/TUFFC.2005.1406545); “Part II: Design and performance for medical imaging applications,” 52(2):192 (doi:10.1109/TUFFC.2005.1406546); “Part III: High frame rate imaging,” 52(2):208 (doi:10.1109/TUFFC.2005.1406547), IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 2005.
R. Y. Chiao, X. Hao, “Coded excitation for diagnostic ultrasound: a system developer's perspective,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 52(2):160–170, 2005.
N. Levanon, E. Mozeson, Radar Signals, Wiley, 2004. Ambiguity functions, phase codes, mismatched filters, NLFM.
S. M. Kay, Fundamentals of Statistical Signal Processing, Vol. II: Detection Theory, Prentice Hall, 1998. Matched filters and generalized matched filters in colored noise.
U.S. FDA, Marketing Clearance of Diagnostic Ultrasound Systems and Transducers, guidance for industry (2019 and later revisions). Track 3 limits and acoustic-output reporting.
IEC 60601-2-37 (basic safety and essential performance of ultrasonic diagnostic equipment) and IEC 62359 (test methods for the determination of thermal and mechanical indices).

← 09 Digital filter design 11 Correlation & delay →