Analytic signal & IQ demodulation
Most ULM localizers, Doppler estimators and phase-based sound-speed methods you will read about run on complex baseband (IQ) data rather than RF (coherence factors are computed on either). This module rebuilds the analytic signal and the IQ pipeline from the spectrum up, and shows exactly what changes in your DAS beamformer when its input is IQ.
You will be able to
- Explain why a real signal's negative frequencies are redundant and build the analytic signal with an FFT, including its edge failures
- Design an IQ chain (mix, low-pass, decimate) and pick \(f_d\), cutoff and decimation factor from the signal bandwidth
- Convert a phase shift to a delay and displacement, and state the wrap limit
- Beamform IQ data correctly (envelope interpolation plus phase rotation) and predict the loss if the rotation is omitted
- Log-compress an envelope and estimate the attenuation-induced center-frequency downshift from IQ
Assumes
- Module 01 · Complex numbers and phasors: phasors, I and Q, demodulation as stopping the rotation, the analytic-signal preview
- Module 03 · Convolution, LTI systems and filters: convolution theorem, low-pass filtering, group delay
- Module 04 · Sampling, Fourier & LTI systems: sampling theorem, aliasing, bandpass sampling, record length and resolution
- DAS beamforming with dynamic receive focusing, plane-wave transmit
08.1Real signals carry their spectrum twice
Take any real signal \(x(t)\) and its Fourier transform \(X(f) = \int x(t)\, e^{-j2\pi f t}\, dt\). Conjugating both sides and using \(x^* = x\) gives the Hermitian symmetry that drives this whole module:
- \(x(t)\)
- real signal (e.g. an RF line)
- \(X(f)\)
- its Fourier transform, a function of frequency \(f\) (MHz)
- \(X^{*}\)
- complex conjugate of \(X\)
- \(\lvert X\rvert\)
- magnitude spectrum
- \(\angle X\)
- phase spectrum, rad
- \(\Longleftrightarrow\)
- if and only if
The elementary case is a cosine. Euler's formula splits it into two counter-rotating phasors of half amplitude:
- \(f_0\)
- carrier (center) frequency, MHz
- \(\varphi\)
- phase at \(t = 0\), rad
- \(t\)
- time, µs
- \(j\)
- imaginary unit, \(j^{2}=-1\)
- \(e^{\pm j(\cdot)}\)
- complex exponentials: counter-rotating phasors at \(\pm f_0\)
A phasor tells you amplitude and phase at one instant. A real sample does not: the value \(\cos\theta = 0.5\) is consistent with \(\theta = +60°\) and \(\theta = -60°\). This sign ambiguity is exactly the redundancy in the spectrum. Discarding the negative-frequency phasor resolves it.
Phasors in the complex plane
Exact formulaThink of the real RF trace as the shadow of a rotating phasor on the real axis. The shadow alone does not reveal which way the phasor turns or where it is on the circle. The analytic signal recovers the phasor from its shadow; IQ demodulation then stops it from spinning at \(f_0\) so only the slow envelope and phase changes remain.
08.2The analytic signal and the Hilbert transform
Build the signal that keeps only the positive half of the spectrum, doubled so that its real part is still \(x(t)\):
- \(X(f)\)
- spectrum of the real signal \(x(t)\), as in (08.1)
- \(X_a(f)\)
- spectrum of the analytic signal (negative frequencies removed)
- \(u(f)\)
- unit step: 1 for \(f>0\), 0 for \(f<0\)
- \(\operatorname{sgn}(f)\)
- sign function: +1 for \(f>0\), −1 for \(f<0\)
- \(x_a(t)\)
- analytic signal (complex), same units as \(x\)
- \(\mathcal{H}\{x\}\)
- Hilbert transform of \(x\): the quadrature component
- \(H(f)\)
- Hilbert transfer function: unit magnitude, −90° for \(f > 0\), +90° for \(f < 0\)
- \(\mathrm{p.v.}\)
- Cauchy principal value of the integral at \(\tau = t\)
- \(x(t)\)
- real signal, units of \(x\)
- \(\leftrightarrow\)
- Fourier-transform pair
- \(\operatorname{sgn}\)
- sign function, (08.3)
- \(j\)
- imaginary unit
- \(\tau\)
- integration variable (time), µs
- \(f\)
- frequency, MHz
- \(\int\)
- integral over all time
Check on the cosine: \(\mathcal{H}\{\cos(2\pi f_0 t)\} = \sin(2\pi f_0 t)\), so \(x_a = \cos + j\sin = e^{j2\pi f_0 t}\), the single phasor of Demo 08.1. Three quantities follow from the polar form \(x_a = A(t)\,e^{j\varphi(t)}\):
- \(A(t)\)
- envelope, same units as \(x\)
- \(\varphi(t)\)
- instantaneous phase, rad
- \(f_i(t)\)
- instantaneous frequency (derivative of the unwrapped phase), MHz
- \(x_a(t)\)
- analytic signal, (08.4)
- \(x(t)\)
- real signal
- \(\lvert\cdot\rvert\)
- magnitude of a complex number
- \(\arg\)
- phase angle of a complex number
- \(t\)
- time, µs
- \(d\varphi/dt\)
- time derivative of the unwrapped phase
Is the envelope the "true" envelope? For \(x(t) = a(t)\cos(2\pi f_0 t)\) the analytic signal equals \(a(t)\,e^{j2\pi f_0 t}\) only if the spectrum of \(a(t)\) lies inside \(\lvert f\rvert < f_0\) (Bedrosian's product theorem). A Gaussian envelope has infinite support in frequency, so the identity is approximate: the error is governed by how much of the positive lobe spills below 0 Hz, a fraction \(\Phi(-f_0/\sigma_f)\) of the Gaussian. At 5 MHz with 60 % bandwidth, \(\sigma_f = 1.27\) MHz and \(f_0/\sigma_f = 3.9\), so the spill is \(4\times10^{-5}\): negligible. At 2 MHz with 120 % bandwidth (\(f_0/\sigma_f = 1.96\)) it is 2.5 %, and the envelope departs from the Gaussian by 2.5 % rms (MATLAB lab, m08_analytic_signal.m).
Mean frequency identity. Weighting the instantaneous frequency by the instantaneous power gives the spectral centroid of the one-sided power spectrum:
- \(f_i(t),\ A(t)\)
- instantaneous frequency and envelope, (08.5)
- \(X_a(f)\)
- analytic spectrum, (08.3)
- \(f\)
- frequency, MHz
- \(t\)
- time, µs
- \(\int_0^\infty\)
- integral over positive frequencies only
- \(\lvert\cdot\rvert^{2}\)
- power (squared magnitude)
Discrete implementation (FFT method)
For an \(N\)-point record (\(N\) even), take the DFT, multiply by a mask and invert:
- \(h[k]\)
- frequency-domain weights that keep positive frequencies (×2) and drop negative ones
- \(k\)
- DFT bin index
- \(N\)
- DFT length (even), samples
- \(X[k]\)
- DFT of the real sequence \(x[n]\)
- \(x[n]\)
- real sampled signal
- \(n\)
- sample index
- \(x_a[n]\)
- analytic signal samples (complex)
- \(\mathrm{IDFT}\)
- inverse discrete Fourier transform
- \(k=0,\ N/2\)
- DC and Nyquist bins (kept once)
The mask is multiplication in the DFT domain, hence circular convolution in time with the discrete Hilbert kernel \(h_{\mathcal{H}}[n] = 2/(\pi n)\) for odd \(n\), 0 for even \(n\). That kernel decays only as \(1/n\), so an echo truncated by the record edge rings across the record, and energy at one end leaks into the other end.
Hilbert transform, envelope and instantaneous frequency of an RF echo
SimulationFFT-mask analytic signal (08.7) of a synthetic 40 MHz record, 1024 samples, computed in the page (JS); speckle line = sum of 160 seeded random Gaussian echoes.Instantaneous frequency is undefined at envelope nulls. In speckle, the phase jumps by up to π where \(\lvert x_a\rvert\) passes near zero, so \(f_i\) spikes far outside the transducer band (even negative). Never average \(f_i\) unweighted; weight by \(A^2\) (which gives the spectral centroid) or use the lag-one autocorrelation estimator of Section 08.6. The same applies to phase-based displacement: phase is unreliable where the echo is weak.
Applying hilbert() to a beamformed frame along the wrong dimension (lateral instead of axial) or to a record that cuts an echo at the edge corrupts the first and last few hundred nanoseconds. Pad with zeros (at least the pulse length) or taper the record before the FFT, and discard the edges.
The B-mode envelope is \(\lvert x_a\rvert\); ULM localizes bubbles on \(\lvert \mathrm{IQ}\rvert\), which is the same quantity (Section 08.3). The axial Hilbert transform is also the reason a beamformed RF image and a beamformed IQ image give the same envelope: both are the same analytic signal, one built after beamforming and one before (Section 08.5).
08.3IQ demodulation: mix, low-pass, decimate
A pulse-echo RF signal is narrowband about \(f_0\). Write it with its complex envelope \(b(t) = A(t)\,e^{j\varphi(t)}\):
- \(b(t)\)
- complex envelope (baseband echo), same units as \(x\)
- \(A(t),\ \varphi(t)\)
- its magnitude and phase (slow compared with \(1/f_0\))
- \(x(t)\)
- real RF signal
- \(t\)
- time, µs
- \(\mathrm{Re}\{\cdot\}\)
- real part
- \(j\)
- imaginary unit
- \(f_0\)
- carrier frequency, MHz, as in (08.2)
Multiplying by \(e^{-j2\pi f_d t}\) (the demodulation frequency \(f_d\) is a design choice, nominally \(f_0\)) shifts the whole spectrum left by \(f_d\):
- \(x(t),\ b(t),\ f_0\)
- RF signal, complex envelope and carrier, (08.8)
- \(f_d\)
- mixing (demodulation) frequency, MHz
- \((\cdot)^{*}\)
- complex conjugate
- \(j\)
- imaginary unit
- \(t\)
- time, µs
- \(f_d\)
- demodulation (mixing) frequency, MHz; stored with the data
- \(\Delta\)
- offset of the echo center from \(f_d\), MHz
- \(\mathrm{LPF}\)
- low-pass filter that passes the baseband lobe and rejects the image at \(-(f_0+f_d)\)
- \(\mathrm{IQ}\)
- \(I + jQ\), with \(I = \mathrm{Re}\{\mathrm{IQ}\}\), \(Q = \mathrm{Im}\{\mathrm{IQ}\}\)
- \(x(t),\ b(t),\ f_0\)
- RF signal, complex envelope and carrier, (08.8)
- \(x_a(t)\)
- analytic signal, (08.4)
- \(t\)
- time, µs
- \(j\)
- imaginary unit
- \(2\)
- compensates the factor ½ in (08.9)
Decimation limit. After the LPF the complex signal occupies roughly \([\Delta - B/2,\ \Delta + B/2]\), where \(B\) is its (two-sided, complex) bandwidth. A complex sample stream at rate \(f_{s,\mathrm{IQ}}\) represents frequencies in \([-f_{s,\mathrm{IQ}}/2,\ +f_{s,\mathrm{IQ}}/2)\) without ambiguity, because there is no mirror image to collide with. So the requirement is
- \(f_s\)
- RF sampling rate, MHz
- \(M\)
- decimation factor (integer)
- \(f_{s,\mathrm{IQ}}\)
- complex (IQ) sample rate, MHz
- \(B\)
- full width of the echo spectrum at the level you must preserve (−20 to −40 dB, not −6 dB), MHz
- \(f_{\max}\)
- highest RF frequency at that level, MHz
- \(\lvert\Delta\rvert\)
- magnitude of the offset in (08.10), MHz
- \(\gt\)
- requirement: rate must exceed the right-hand side to avoid aliasing
LPF cutoff tradeoff. The cutoff \(f_c\) must pass \(\lvert\Delta\rvert + B/2\) (or the envelope is smoothed and the phase biased) and must reject the image at \(-(f_0+f_d)\) and the noise above \(f_{s,\mathrm{IQ}}/2\) (or it aliases into the band after decimation). With a 49-tap Hamming FIR at 40 MHz the transition band is \(\approx 3.3 f_s/L \approx 2.7\) MHz wide, which is why \(f_c\) needs margin on both sides.
The IQ pipeline, stage by stage
SimulationGaussian echo sampled at 40 MHz, complex mixer, 49-tap Hamming-windowed-sinc FIR and decimation, all computed in the page (JS); reference envelope from the FFT Hilbert mask.Aliasing at decimation does not change the samples you keep. It corrupts everything you compute between them: fractional-delay interpolation in the beamformer, upsampling for display, axial phase-based delay estimation and any spectral estimate. A decimated IQ set can look fine in a sample-by-sample comparison and still beamform badly.
IQ exports and research pipelines usually fix \(f_d\) per acquisition, commonly at or near the nominal transmit frequency rather than the received center frequency (Section 08.7 covers concrete pipelines). Because attenuation lowers the received center frequency with depth (Section 08.6), real data always has some depth-dependent \(\Delta\). The envelope is unaffected; anything converting phase to time must account for it (Section 08.4).
08.4Phase carries delay
Delay the echo by \(\tau\). With \(f_d = f_0\), the IQ signal of \(x(t - \tau)\) is
- \(\mathrm{IQ}(t)\)
- baseband IQ signal, (08.10)
- \(\mathrm{IQ}_\tau(t)\)
- IQ signal delayed by \(\tau\)
- \(b(t)\)
- complex envelope, (08.8)
- \(\tau\)
- time delay, µs
- \(f_0\)
- carrier frequency, MHz
- \(t\)
- time, µs
- \(j\)
- imaginary unit
- \(\ll\)
- much smaller than
- \(\Delta\varphi\)
- phase of the zero-lag complex correlation of two IQ records, rad, in \((-\pi, \pi]\)
- \(\hat{\tau}\)
- estimated two-way delay of record 2 relative to record 1, s
- \(f_0\)
- conversion frequency; must be the echo's mean frequency, Hz
- \(\hat{d}\)
- axial displacement, m (positive away from the probe)
- \(c\)
- speed of sound, 1540 m/s
- \(\mathrm{IQ}_1,\ \mathrm{IQ}_2\)
- IQ records before and after the delay or displacement
- \(n\)
- fast-time sample index
- \(\sum_n\)
- sum over the samples in the window
- \((\cdot)^{*}\)
- complex conjugate
- \(\arg\)
- phase angle
- \(\Rightarrow\)
- implies
| Quantity (\(f_0\) = 5 MHz, \(c\) = 1540 m/s) | Expression | Value |
|---|---|---|
| Delay per radian | \(1/(2\pi f_0)\) | 31.8 ns |
| Round-trip displacement per radian | \(c/(4\pi f_0) = \lambda/(4\pi)\) | 24.5 µm |
| Phase for a 10 ns delay | \(2\pi f_0 \tau\) | 0.314 rad = 18° |
| Unambiguous delay range | \(\lvert\tau\rvert < 1/(2f_0)\) | ±100 ns |
| Unambiguous axial displacement | \(\lvert d\rvert < \lambda/4\) | ±77 µm |
| Kasai Nyquist velocity, PRF 5 kHz | \(c\cdot\mathrm{PRF}/(4f_0)\) | 0.385 m/s |
Two systematic errors follow directly from the conversion step. First, the phase wraps: a delay of \(1/f_0 - 10\) ns is indistinguishable from −10 ns. Second, \(\hat{\tau}\) uses a frequency, and the right one is the mean frequency of the echo, not the nominal \(f_0\) or \(f_d\). If the true mean frequency is 10 % below the value used, every delay is underestimated by 10 %. Attenuation causes exactly this (Section 08.6).
Phase-based delay estimate and its wrap
SimulationTwo sampled Gaussian echoes, IQ formed and correlated in the page (JS); estimate from (08.13).Module 11 treats phase-based delay estimators properly: the lag-one autocorrelator of Kasai et al. (slow time, Doppler), its 2-D extension by Loupas et al. (which estimates the mean RF frequency from the fast-time lag and so removes the conversion bias above), and their variance versus cross-correlation. Sound-speed target: computed ultrasound tomography in echo mode (CUTE, Jaeger et al.) measures the phase shift between IQ images beamformed from different steering angles; a wrong assumed \(c\) produces angle-dependent echo delays, so the phase map \(\Delta\varphi(x,z) = -2\pi f_0\, \Delta\tau(x,z)\) is inverted for \(c(x,z)\). Aberration target: a phase screen is a per-element delay error, estimated as the phase of the correlation between neighboring channels' IQ signals and wrapped at \(\lambda/2\) of path difference.
08.5Beamforming on IQ
You know DAS on RF: sample channel \(n\) at the focusing delay \(\tau_n(\mathbf{r})\) and sum. On IQ, the carrier has been removed, and that is what changes. Channel \(n\) receives an echo from a scatterer at round-trip time \(t_n\). Demodulated with \(f_d\):
- \(t_n\)
- round-trip time of the echo on channel \(n\), s
- \(\tau_n(\mathbf{r})\)
- focusing delay of channel \(n\) for image point \(\mathbf{r}\), s
- \(\mathrm{IQ}_n(t)\)
- baseband IQ signal of channel \(n\)
- \(b(t)\)
- complex envelope, (08.8)
- \(f_d,\ \Delta\)
- demodulation frequency and offset, (08.10)
- \(t\)
- time, µs
- \(j\)
- imaginary unit
- \(n\)
- channel (element) index
Sampling at the focusing delay \(\tau_n = t_n\) gives \(b(0)\,e^{-j2\pi f_d \tau_n}\). The envelope is aligned, but each channel carries a different phase. The correct IQ beamformer therefore applies the phase rotation that the RF carrier would have applied automatically:
- \(y_{\mathrm{IQ}}(\mathbf{r})\)
- beamformed complex pixel at \(\mathbf{r} = (x, z)\)
- \(w_n\)
- apodization weight of channel \(n\)
- \(\mathrm{IQ}_n(\tau_n)\)
- channel IQ interpolated at the focusing delay (slow envelope, so linear interpolation at \(f_{s,\mathrm{IQ}}\) is adequate)
- \(e^{+j2\pi f_d\tau_n}\)
- phase rotation by the full delay times the demodulation frequency \(f_d\), not \(f_0\)
- \(\mathbf{r}\)
- image point \((x,z)\), mm
- \(n\)
- channel index
- \(\sum_n\)
- sum over channels
- \(f_d\)
- demodulation frequency, MHz, (08.10)
- \(j\)
- imaginary unit
This is not a narrowband approximation. Since \(\mathrm{IQ}_n(t) = x_{a,n}(t)\,e^{-j2\pi f_d t}\) holds identically, \(\mathrm{IQ}_n(\tau)\,e^{+j2\pi f_d \tau} = x_{a,n}(\tau)\) for every \(\tau\) and every frequency in the band. The only errors are those of the IQ interpolator and the demodulation LPF. In the MATLAB lab (32 elements, 5 MHz, 49-tap LPF) they cost 0.4 dB of peak with linear interpolation at 10 MHz and 0.1 dB at 20 MHz. A beamformer that applies only the phase rotation and skips the envelope shift (phase-shift beamforming) is the narrowband approximation; it fails once the delay spread across the aperture approaches the pulse length.
Omitting the rotation sums \(b(0)\,e^{-j2\pi f_d \tau_n}\) across the aperture. With \(\tau_n \approx (2z + x_n^2/(2z))/c\) (Fresnel approximation, plane-wave transmit), the phase across the aperture is quadratic with peak excursion \(\psi = \pi x_{\max}^2/(\lambda z)\):
- \(N\)
- number of active elements
- \(x_{\max}\)
- lateral position of the outermost element from the aperture center, mm
- \(\lambda\)
- wavelength \(c/f_d\), mm
- \(z\)
- target depth, mm
- \(s\)
- normalized aperture coordinate \(x/x_{\max}\)
- \(G_{\text{no rot}}\)
- coherent gain when the \(e^{j2\pi f_d\tau_n}\) rotation of (08.15) is omitted (0 to 1, dimensionless)
- \(f_d\)
- demodulation frequency, MHz
- \(\tau_n\)
- focusing delay of channel \(n\), µs
- \(n\)
- channel index
- \(\sum_n\)
- sum over channels
- \(\psi\)
- quadratic phase across the aperture, rad
- \(j\)
- imaginary unit
- \(\int_0^1\)
- continuous approximation of the sum over the aperture
- \(\lvert\cdot\rvert\)
- magnitude
- \(\approx\)
- Fresnel (parabolic delay) approximation
DAS on RF vs IQ, with and without phase rotation
SimulationPoint-target channel data from straight-ray two-way delays (0° plane wave, Gaussian 60 % pulse, no diffraction or element directivity), demodulated, decimated and beamformed with (08.15) in the page (JS).Rotating by \(e^{j2\pi f_0\tau}\) with the nominal transmit \(f_0\) when the data were demodulated at a different \(f_d\) leaves a residual \(e^{j2\pi(f_d - f_0)\tau_n}\) across the aperture. Use the demodulation frequency stored with the data. The same rule applies to any IQ re-delay: synthetic aperture, retrospective transmit focusing, coherent compounding of IQ frames with different steering delays, and motion compensation.
Most ULM, functional ultrasound and ultrafast Doppler pipelines beamform IQ, because the compounded IQ frames feed the SVD clutter filter (Module 06) and the slow-time phase is the Doppler signal. A missing or wrong rotation shows up as a depth-dependent loss of SNR and a broadened PSF, which directly worsens localization precision in ULM. In the differentiable-beamforming sound-speed methods, the phase term \(e^{j2\pi f_d \tau_n(c)}\) is the part of the forward model most sensitive to \(c\): a delay error of \(1/(2\pi f_d)\) is one radian of incoherence.
08.6Envelope, log compression and the attenuation downshift
Display. The B-mode value is the envelope \(\lvert\mathrm{IQ}\rvert\), log-compressed into a fixed display window of dynamic range DR:
- \(A_{\mathrm{ref}}\)
- reference amplitude mapped to white (0 dB)
- \(\mathrm{DR}\)
- displayed dynamic range, dB
- \(I_{\mathrm{dB}}\)
- displayed intensity, dB (amplitude, hence 20 log)
- \(\mathrm{IQ}\)
- beamformed complex pixel, (08.15)
- \(\lvert\cdot\rvert\)
- envelope (magnitude)
- \(\log_{10}\)
- base-10 logarithm
- \(\mathrm{gray}\)
- display gray level, 0 (black) to 1 (white)
- \(\mathrm{clip}(\cdot,0,1)\)
- limits the value to the range 0 to 1
Rayleigh teaser. For fully developed speckle (many random scatterers per resolution cell), IQ is circular complex Gaussian, so \(\lvert\mathrm{IQ}\rvert\) is Rayleigh distributed. Its mean-to-standard-deviation ratio is \(\sqrt{\pi/(4-\pi)} = 1.91\) independent of scatterer strength, and the standard deviation of \(20\log_{10}\lvert\mathrm{IQ}\rvert\) is 5.57 dB. Speckle statistics get their own module (Module 26).
Frequency downshift. Let the transmitted amplitude spectrum be Gaussian with center \(f_0\) and standard deviation \(\sigma_f\) (\(\sigma_f = \mathrm{BW}_{-6\,\mathrm{dB}}/2.355\), with \(\mathrm{BW}_{-6\,\mathrm{dB}}\) the −6 dB width of the amplitude spectrum). The amplitude attenuation coefficient is linear in frequency, \(\alpha(f) = \alpha_0 f\) with \(\alpha_0\) in Np/cm/MHz, and the echo from one-way depth \(z\) travels \(2z\), so the amplitude spectrum is multiplied by \(e^{-2\alpha_0 f z}\). Complete the square in the exponent:
- \(f\)
- frequency, MHz
- \(f_0\)
- centre frequency of the transmitted pulse, MHz
- \(\sigma_f\)
- standard deviation of the Gaussian amplitude spectrum, MHz
- \(\alpha_0\)
- attenuation slope, Np/(cm·MHz), (08.19)
- \(z\)
- one-way depth, cm
- \(\mathrm{const}(z)\)
- terms that depend on depth but not on frequency
- \(f_c(z)\)
- center frequency of the echo from depth \(z\), MHz
- \(\alpha_0\)
- amplitude attenuation slope, Np/(cm·MHz) (divide dB values by 8.686)
- \(\sigma_f\)
- standard deviation of the amplitude spectrum, MHz
- \(\sigma_p\)
- \(\sigma_f/\sqrt{2}\): standard deviation of the power spectrum \(\lvert S\rvert^2\), MHz
- \(z\)
- one-way depth, cm
- \(f_0\)
- centre frequency of the transmitted pulse, MHz
- \(\Longleftrightarrow\)
- equivalent form in terms of the power-spectrum width
Worked number. \(f_0\) = 5 MHz, 60 % bandwidth → \(\sigma_f = 3/2.355 = 1.274\) MHz, \(\sigma_f^2 = 1.623\) MHz². Liver-like \(\alpha_0\) = 0.5 dB/cm/MHz = 0.0576 Np/cm/MHz. Slope \(-2\alpha_0\sigma_f^2 = -0.187\) MHz/cm; at \(z\) = 4 cm, \(\Delta f = -0.75\) MHz (15 %). A simulation that applies \(e^{-2\alpha_0 f z}\) to the Gaussian and takes the power-spectrum centroid gives −0.747 MHz. Inverting the measured slope gives the spectral-shift attenuation estimate \(\hat{\alpha}_0 = -(df_c/dz)/(2\sigma_f^2)\).
Estimating \(f_c\) from IQ. Apply the mean-frequency identity (08.6) in discrete form. With IQ demodulated at \(f_d\) and sampled at \(f_{s,\mathrm{IQ}}\), the lag-one autocorrelation in fast time over a depth window \(W\) gives
- \(W\)
- depth window (fast-time samples), typically a few mm
- \(R(1)\)
- lag-one fast-time autocorrelation, complex
- \(\hat{f}_c\)
- local mean (center) frequency estimate, MHz
- \(\mathrm{IQ}[k]\)
- fast-time IQ samples along depth, (08.10)
- \(k\)
- fast-time sample index
- \(\sum_{k\in W}\)
- sum over the depth window \(W\)
- \((\cdot)^{*}\)
- complex conjugate
- \(f_d\)
- demodulation frequency, MHz
- \(f_{s,\mathrm{IQ}}\)
- IQ sample rate, MHz, (08.11)
- \(\arg\)
- phase angle, rad
B-mode log compression, Rayleigh envelope and depth-dependent downshift
SimulationCircular Gaussian scatterer field at 20 MHz IQ rate convolved with a depth-dependent Gaussian baseband pulse (spectrum × \(e^{-2\alpha_0 f z}\)), 64 lines, computed in the page (JS); 1D convolution model, no diffraction.This is the spectral-shift branch of QUS attenuation imaging. In practice the measured shift is contaminated by diffraction (the focal gain is frequency dependent), by the depth-dependent \(f_d/f_0\) mismatch of the IQ export, and by non-Gaussian transducer spectra; reference-phantom methods divide these out. The same downshift biases phase-to-delay conversion (Section 08.4) and the \(f_d\) used for IQ beamforming phase, so deep phase-based sound-speed and displacement estimates should use a locally estimated \(\hat{f}_c\).
Attenuation is a low-pass filter whose slope grows with depth. Multiply a Gaussian by a decaying exponential and you get a Gaussian shifted toward the decaying side; nothing else changes. That is why one number (the slope) summarizes the effect, and why it scales with \(\sigma_f^2\).
08.7IQ in practice: research pipelines and their choices
Three decisions define any IQ pipeline: where the complex signal is formed (per channel or after beamforming), at what rate it is stored, and which demodulation frequency is used. Each sets what the downstream estimators may assume.
| Pipeline | Where IQ is formed | Beamforms | Why this way | What it costs |
|---|---|---|---|---|
| Conventional (hardware) scanner | After beamforming, per image line | RF, with fine delays from oversampling or interpolation filters | One demodulator per line instead of per channel; depth-dependent filtering is easy on one line | RF delay accuracy needs a high effective rate (below) |
| Software beamformer on a research platform | Per channel, before beamforming | IQ, with interpolation + phase rotation | Low channel data rate on the bus and GPU; exact fine delay comes from the rotation, not the sample grid; Doppler, SVD and ULM consume IQ anyway | \(f_d\) must be stored and used; one complex multiply per channel sample |
| Quadrature (\(4f_d\)) sampling | In the sampling itself: \(f_s = 4f_d\) | Either | The mixer \(e^{-j\pi n/2}\) is (1, −j, −1, +j): sign flips only, no multiplies | I and Q come from samples \(1/(4f_d)\) apart; needs interpolation for broadband pulses (Demo 08.7) |
Four samples per wavelength. Verasonics-style systems tie the RF sample rate to a demodulation frequency, commonly 4 samples per wavelength (\(f_s = 4f_d\)), and form IQ in software. With \(f_s = 4f_d\) the sampled carrier is \(\cos(\pi n/2 + \varphi)\) and the samples cycle through \(+A\cos\varphi\), \(-A\sin\varphi\), \(-A\cos\varphi\), \(+A\sin\varphi\): even samples carry \(\pm I\), odd samples carry \(\mp Q\). Four samples per wavelength is real sampling at \(4f_d\), enough for any echo whose spectrum stays inside \((0, 2f_d)\), a total width of \(2f_d\). For a Gaussian pulse centered at \(f_d\) whose −40 dB edges must fit, that allows a −6 dB fractional bandwidth up to about 78 %. Reduced modes keep fewer samples per wavelength to save data and accept less bandwidth.
- \(x[n]\)
- real samples taken at \(f_s=4f_0\)
- \(n\)
- sample index
- \(k\)
- block index (groups of four samples)
- \(A,\ \varphi\)
- amplitude and phase of the carrier
- \(f_0\)
- carrier frequency (quarter of the sampling rate)
- \(I,\ Q\)
- in-phase and quadrature components
- \(j\)
- imaginary unit
- \(\Rightarrow\)
- implies
IQ rate versus bandwidth. Two reals per complex sample means IQ at \(f_s/4\) carries the same number of reals as RF at \(f_s/2\). The saving is relative to oversampled RF, not to RF at its Nyquist rate. The real advantage of channel IQ is that the carrier phase is restored analytically by \(e^{j2\pi f_d \tau}\), so the interpolator only has to follow the slow envelope. RF DAS must interpolate the carrier itself: with nearest-sample delays the phase error is uniform in \(\pm\pi f_0/f_s\) and the coherent gain is \(\sin(a)/a\) with \(a = \pi f_0/f_s\), −0.22 dB at 40 MHz and 5 MHz, plus quantization sidelobes. Hence RF beamformers oversample or interpolate.
Choosing \(f_d\) when the spectrum downshifts. Attenuation moves the echo center from \(f_0\) toward \(f_0 - 2\alpha_0\sigma_f^2 z\) (08.19), 0.75 MHz at 4 cm in liver at 5 MHz, 60 %. A fixed \(f_d = f_0\) then leaves a depth-growing \(\Delta\), which must fit under the LPF and inside \(f_{s,\mathrm{IQ}} > B + 2\lvert\Delta\rvert\) (08.11). Options: put \(f_d\) at the center frequency expected at mid-depth of interest; or let \(f_d\) track depth. If \(f_d\) varies, the demodulator phase is \(\theta_d(t) = 2\pi\int_0^t f_d(t')\,dt'\), and the IQ beamformer must rotate by \(e^{j\theta_d(\tau)}\), which equals \(e^{j2\pi f_d \tau}\) only for constant \(f_d\). In every case: rotate with what the demodulator removed, convert phase to delay with the local echo frequency \(\hat{f}_c\).
RF or IQ beamforming? Beamform RF when you need the full real spectrum after beamforming (harmonic or spectral QUS processing that wants the band untouched by a demodulation LPF) or when the hardware already does it. Beamform IQ when data rate and GPU memory dominate (ultrafast plane-wave, ULM, functional ultrasound), when the outputs feed slow-time phase processing, or when you will re-beamform the same channel data many times (sound-speed sweeps, autofocus): the rotation makes re-delaying exact at a low sample rate.
IQ samples of one short echo at fs = 4f0
SimulationGaussian echo sampled at 20 MHz, IQ formed by sample pairing or by the (1, −j, −1, +j) mixer and [¼ ½ ¼] filter, computed in the page (JS) and compared with the analytic baseband.IQ data from different systems differ in sign convention (\(e^{-j2\pi f_d t}\) versus \(e^{+j2\pi f_d t}\), i.e. \(I - jQ\) versus \(I + jQ\)), in the factor 2 after the LPF, and in whether \(f_d\) is stored at all. Before any phase-based processing, check with one synthetic echo: delay it by a known \(\tau\) and confirm the measured phase is \(-2\pi f_c \tau\) with your sign, and that the Doppler of a target moving toward the probe comes out positive.
08.8Visual summary
Eleven pictures that carry the module, each computed from its formula: a question you should answer in five seconds, the answer with its number, and the figure to keep in mind. The slider sets the envelope width of the transform-pair card: a longer pulse always means a narrower lobe.
08.9How experts think about it
An expert sees every RF echo as a slow complex envelope riding a known carrier. The envelope (magnitude) is a property of the echo; the phase is a property of the echo and of the bookkeeping (which \(f_d\), which sign, which time origin). Every IQ bug is a bookkeeping bug, so the reflex is to ask which frequency was removed, which frequency the echo actually has, and whether the two are being confused.
\(\mathrm{IQ}(t) = x_a(t)\,e^{-j2\pi f_d t}\) with an ideal LPF (08.10), and \(\lvert e^{j\theta}\rvert = 1\). So B-mode, ULM detection and envelope statistics are identical for any \(f_d\) that keeps the band inside the LPF and the IQ Nyquist range. The phase is shifted by the deterministic ramp \(2\pi f_d t\), which must be known to reconstruct RF, re-beamform or compare data demodulated differently.
A complex stream has no mirror image, so its spectrum may occupy the full \([-f_s/2, f_s/2)\). Real sampling must keep the positive band away from its mirror, hence \(f_s > 2f_{\max}\). For a 5 MHz, 80 % pulse: IQ at 8 to 10 MHz complex versus RF ≥ 17.3 MHz at the −20 dB edge (8.65 MHz), 30 to 40 MHz in practice.
After demodulation, channel \(n\) sampled at its focusing delay carries \(b(0)\,e^{-j2\pi f_d \tau_n}\). Without rotation the aperture sum is a Fresnel integral \(\bigl\lvert\int_0^1 e^{-j\psi s^2} ds\bigr\rvert\) (08.16). 32 elements, 0.3 mm pitch, 5 MHz, \(z\) = 20 mm: \(\psi\) = 11.0 rad, Fresnel −12.5 dB, exact sum −13.4 dB, simulated −13.3 dB. The loss oscillates with depth (−14 to −6.4 dB from 10 to 60 mm) and the PSF smears laterally over millimetres.
\(\mathrm{IQ}_n(\tau)\,e^{j2\pi f_d \tau} = x_{a,n}(\tau)\) holds identically for every frequency in the band, so rotated IQ DAS equals the analytic signal of RF DAS. Residual errors come from interpolating the baseband envelope and from LPF droop. In m08_iq_beamforming.m (49-tap LPF, \(f_c\) = 4 MHz) the peak is 0.4 dB below RF DAS with linear interpolation at 10 MHz and 0.1 dB below at 20 MHz.
A delay \(\tau\) rotates IQ by \(-2\pi f_0\tau\). Round trip: \(\tau = 2d/c\), so \(d\) per radian \(= c/(4\pi f_0) = \lambda/(4\pi)\) = 0.308 mm/12.57 = 24.5 µm. One degree is 0.43 µm. This is why phase tracking resolves sub-micrometre motion while the envelope shift of 10 ns (0.4 samples at 40 MHz) is invisible in the magnitude.
arg returns \((-\pi, \pi]\), so \(\lvert 2\pi f_0\tau\rvert < \pi\), i.e. \(\lvert\tau\rvert < 1/(2f_0)\), and round trip halves it again: \(\lvert d\rvert < \lambda/4\). Between frames at PRF 5 kHz that is 77 µm × 5000 s⁻¹ = 0.385 m/s, the Kasai Nyquist velocity at 5 MHz.
The beamformer's rotation must undo exactly what the demodulator removed (\(f_d\)). The phase a delay produces is \(-2\pi f_{\mathrm{echo}}\tau\), so converting phase to time with \(f_d\) scales every delay by \(f_{\mathrm{echo}}/f_d\). A mismatch also leaves a residual rotation \(e^{j2\pi\Delta t}\) along fast time that biases any estimator comparing IQ at different fast times. Kasai velocity, CUTE-type sound-speed maps and aberration phase screens are all biased unless the conversion uses a local \(\hat{f}_c\) (e.g. Loupas' 2-D autocorrelator).
\(\exp\bigl(-(f-f_0)^2/2\sigma_f^2\bigr)\) \(\cdot\exp(-2\alpha_0 f z)\) is a Gaussian of the same \(\sigma_f\) centered at \(f_0 - 2\alpha_0\sigma_f^2 z\) (complete the square, (08.18)). 5 MHz, 60 %: \(\sigma_f\) = 1.274 MHz; 0.5 dB/cm/MHz = 0.0576 Np/cm/MHz; slope 0.187 MHz/cm. m08_attenuation_downshift.m matches to \(2\times10^{-4}\) MHz and recovers \(\alpha_0\) = 0.49 dB/cm/MHz from speckle.
\(2\alpha_0\sigma_f^2 z = 4\alpha_0\sigma_p^2 z\): the same physics, with \(\sigma\) of the amplitude or of the power spectrum. Mixing the two conventions is a factor of 2; using dB where Np is required is a factor of 8.686. Units check: Np/(cm·MHz) × MHz² × cm = MHz.
The envelope is an amplitude, so dB = \(20\log_{10}\); \(10\log_{10}\) applies to \(\lvert\mathrm{IQ}\rvert^2\). Echoes span 60 to 80 dB, so the display window clips. Log compression turns multiplicative Rayleigh speckle into an additive term of fixed std 5.57 dB, independent of echogenicity; mean/std of the linear envelope is 1.91.
The mask (08.7) multiplies DFT bins, which is circular convolution with a kernel that decays only as \(1/n\). An echo cut by the record start rings across the record and wraps: in m08_analytic_signal.m a truncated echo leaves \(\lvert x_a\rvert\) = 0.095 (−20 dB) in the last 0.5 µs of a 25.6 µs record; 2× zero padding removes it.
\(f_i = (1/2\pi)\,d\varphi/dt\), and phase noise scales as noise/\(A\). Near envelope nulls (speckle) the phase jumps by up to \(\pi\) and \(f_i\) spikes outside the band. Weight by \(A^2\): the \(A^2\)-weighted mean \(f_i\) equals the spectral centroid (08.6). With 5 % noise, the weighted mean in the echo window is 5.03 MHz, the unweighted 5.09 MHz.
\(\cos(\pi n/2 + \varphi)\) cycles through \(\cos\varphi,\ -\sin\varphi,\ -\cos\varphi,\ \sin\varphi\) (08.21). Pairing neighbors gives I and Q 50 ns apart at 5 MHz, half the \(2f_0\) output period. For a 60 % pulse that costs 18 % rms complex error on average; averaging the two Q neighbors ([¼ ½ ¼] after mixing) brings it to 4 %.
\(a(t)\cos(2\pi f_0 t)\) has analytic signal \(a(t)\,e^{j2\pi f_0 t}\) iff the spectrum of \(a\) lies inside \(\lvert f\rvert < f_0\). A Gaussian lobe spills a fraction \(\Phi(-f_0/\sigma_f)\) below 0 Hz: \(4\times10^{-5}\) at 5 MHz, 60 %; 2.5 % at 2 MHz, 120 %, where the envelope differs from the Gaussian by 2.5 % rms.
Decimation keeps samples unchanged, so a sample-by-sample comparison always passes. Aliasing appears in anything computed between samples: beamforming interpolation, upsampling, axial phase-delay estimation, spectra. For a 5 MHz, 60 % pulse the band-limited reconstruction error stays at the 0.6 % LPF floor down to 8 MHz IQ and grows to 3.9 % at 5 MHz and 12 % at 4.4 MHz.
\(R(1) = \sum \mathrm{IQ}^*[k]\,\mathrm{IQ}[k+1]\) estimates \(\int P(f)\, e^{j2\pi f/f_{s,\mathrm{IQ}}}\, df\). For a Gaussian PSD that is a real Gaussian times \(e^{j2\pi(f_c - f_d)/f_{s,\mathrm{IQ}}}\), so the argument is exact. In fast time it gives the local \(\hat{f}_c\) (spectral-shift QUS, phase-to-delay conversion); in slow time it is Kasai's Doppler estimator.
Rounding delays to the sample grid gives a delay error uniform in \(\pm 1/(2f_s)\), hence a carrier phase error uniform in \(\pm\pi f_0/f_s\). The mean phasor has magnitude \(\sin(a)/a\): \(a = \pi/8\) gives 0.974 (−0.22 dB). The periodic error pattern across the aperture also produces quantization sidelobes. IQ DAS gets the carrier phase analytically, which is why it tolerates 10 MHz sampling.
08.10Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Hermitian symmetry | \(X(-f) = X^*(f)\) for real \(x\) | Why one-sided (IQ) storage loses nothing |
| Hilbert transform | \(H(f) = -j\operatorname{sgn}(f)\); \(\mathcal{H}\{\cos\} = \sin\); \(\mathcal{H}\{\mathcal{H}\{x\}\} = -x\) | Quadrature component of RF |
| Analytic signal | \(x_a = x + j\mathcal{H}\{x\}\); \(X_a = 2u(f)X\) | Envelope, phase, mean frequency |
| FFT mask | \(h = [1, 2 \dots 2, 1, 0 \dots 0]\) (DC, Nyquist ×1) | hilbert(); pad or taper, it is circular |
| Envelope / phase / \(f_i\) | \(\lvert x_a\rvert\), \(\arg x_a\), \(\frac{1}{2\pi}\frac{d\varphi}{dt}\) | B-mode, ULM detection, displacement |
| Bedrosian condition | \(a(t)\cos\omega_0 t \mapsto a\,e^{j\omega_0 t}\) iff \(\operatorname{supp} A(f) \subset \lvert f\rvert < f_0\) | Envelope error for very broadband, low-\(f_0\) pulses |
| Mean-frequency identity | \(\langle f_i\rangle_{A^2}\) = centroid of \(\lvert X_a\rvert^2\) | Kasai / centroid estimators |
| IQ demodulation | \(\mathrm{IQ} = 2\,\mathrm{LPF}\{x\,e^{-j2\pi f_d t}\} = x_a\,e^{-j2\pi f_d t}\) | Scanner IQ export |
| Complex sampling limit | \(f_{s,\mathrm{IQ}} > B + 2\lvert f_0 - f_d\rvert\) | \(f_s/4\) decimation is typical |
| Quadrature sampling | \(f_s = 4f_d\): \(x[4k \dots 4k+3] = I, -Q, -I, +Q\) | Multiplier-free IQ; mind the \(1/(4f_d)\) I/Q offset |
| Delay ↔ phase | \(\Delta\varphi = -2\pi f_c\tau\); 1 rad = \(1/(2\pi f_c)\) = \(\lambda/(4\pi)\) axial; wrap at \(\pm 1/(2f_c)\), \(\pm\lambda/4\) | Displacement, Doppler, CUTE, phase screens |
| Kasai Nyquist velocity | \(v_N = c\cdot\mathrm{PRF}/(4f_c)\) | 0.385 m/s at 5 MHz, 5 kHz |
| IQ DAS | \(\sum w_n\,\mathrm{IQ}_n(\tau_n)\,e^{j2\pi f_d\tau_n}\) (variable \(f_d\): \(e^{j\theta_d(\tau_n)}\)) | Most ultrafast / ULM beamformers |
| Loss without rotation | \(\bigl\lvert\frac{1}{N}\sum e^{-j2\pi f_d\tau_n}\bigr\rvert \approx \bigl\lvert\int_0^1 e^{-j\psi s^2}ds\bigr\rvert\), \(\psi = \pi x_{\max}^2/(\lambda z)\) | Diagnosing a broken IQ beamformer |
| RF delay quantization | gain \(\sin(a)/a\), \(a = \pi f_0/f_s\) | Why RF beamformers oversample or interpolate |
| Log compression | \(20\log_{10}(\lvert\mathrm{IQ}\rvert/A_{\mathrm{ref}})\), DR 40–60 dB | B-mode display |
| Rayleigh speckle | mean/std = 1.91; std\((20\log A)\) = 5.57 dB | Fully developed speckle check |
| Attenuation downshift | \(f_c(z) = f_0 - 2\alpha_0\sigma_f^2 z = f_0 - 4\alpha_0\sigma_p^2 z\) (\(\alpha_0\) in Np/cm/MHz, \(z\) one-way) | QUS spectral-shift attenuation |
| Fast-time centroid | \(\hat{f}_c = f_d + \frac{f_{s,\mathrm{IQ}}}{2\pi}\arg\sum \mathrm{IQ}^*[k]\,\mathrm{IQ}[k+1]\) | Local \(f_c\) for phase-to-delay conversion |
08.11Interview and qualifying-exam questions
I1What is \(H(f)\) for the Hilbert transform, what is \(\mathcal{H}\{\cos\}\) and \(\mathcal{H}\{\sin\}\), and what is \(\mathcal{H}\) applied twice? · screen
\(H(f) = -j\operatorname{sgn}(f)\): −90° on positive frequencies, +90° on negative, unit magnitude. \(\mathcal{H}\{\cos\omega_0 t\} = \sin\omega_0 t\), \(\mathcal{H}\{\sin\omega_0 t\} = -\cos\omega_0 t\). Twice: \(H^2 = (-j\operatorname{sgn} f)^2 = -1\) for \(f \ne 0\), so \(\mathcal{H}\{\mathcal{H}\{x\}\} = -x\) (for \(x\) with no DC). The analytic signal \(x + j\mathcal{H}\{x\}\) has spectrum \(2u(f)X(f)\). Follow-up: why is the time-domain kernel \(1/(\pi t)\) a problem in practice? It decays slowly, so FIR Hilbert transformers need many taps and FFT implementations are circular.
I2Why can complex IQ be sampled at a rate set by bandwidth while RF needs twice its highest frequency? Give numbers for a 5 MHz, 80 % probe. · screen
A real signal's spectrum has a mirror at negative frequencies; sampling must keep the positive band from overlapping the mirror's replicas, hence \(f_s > 2f_{\max}\). IQ has only one band, so \(f_{s,\mathrm{IQ}} > B\) (plus \(2\lvert\Delta\rvert\) for a demodulation offset). \(\sigma_f = 4/2.355 = 1.70\) MHz; the −20 dB width is 7.3 MHz, so 8 to 10 MHz complex; RF needs ≥ 17.3 MHz at the −20 dB edge, 30 to 40 MHz in practice. Follow-up: is IQ cheaper in bits? At \(f_s/4\) complex it carries the same number of reals as RF at \(f_s/2\); the gain is against oversampled RF and in interpolation accuracy, not against Nyquist-rate RF.
I3Three ways to get a B-mode envelope from RF. Which are exact? · screen
(1) \(\lvert x_a\rvert\) via Hilbert/FFT: exact by definition (subject to edge effects). (2) IQ demodulation then \(\lvert\mathrm{IQ}\rvert\): equal to (1) for an ideal LPF that passes the band and rejects the image; a real FIR adds sub-percent error (0.6 % in the MATLAB lab). (3) Rectify and low-pass: an approximation; \(\lvert\cos\rvert\) has mean \(2/\pi\) and harmonics at \(2f_0\), so the result depends on the LPF and is scaled by \(2/\pi\). Follow-up: do (1) and (2) differ in phase? Yes, by \(2\pi f_d t\).
I4A colleague's IQ beamformer gives a PSF that is 13 dB weak and smeared over several millimetres at 20 mm, but fine B-mode speckle texture. Diagnose. · onsite
Missing (or wrong-sign) phase rotation. Each channel carries \(e^{-j2\pi f_d\tau_n}\) after envelope alignment; summing without \(e^{+j2\pi f_d\tau_n}\) gives a Fresnel-integral loss, −13.4 dB for 32 elements, 0.3 mm pitch, 5 MHz at 20 mm, oscillating with depth. Speckle still looks like speckle because random-phase sums have the same statistics. Test: beamform a point target at several depths; a loss that oscillates with depth and tracks \(\bigl\lvert\frac{1}{N}\sum e^{-j2\pi f_d\tau_n}\bigr\rvert\) confirms it. A wrong-sign rotation doubles the phase error instead of cancelling it. Follow-up: which \(f\) to use in the rotation? The stored \(f_d\), not the nominal \(f_0\).
I5With 5 MHz IQ, what axial displacement corresponds to 1° of phase, and what is the largest inter-frame motion you can track by phase? · onsite
\(d\) per radian \(= \lambda/(4\pi)\) = 24.5 µm, so 1° = 0.43 µm. Phase wraps at \(\pm\pi\), i.e. \(\pm\lambda/4\) = ±77 µm between frames; at PRF 5 kHz that is 0.385 m/s (Kasai Nyquist \(c\cdot\mathrm{PRF}/(4f_0)\)). Follow-up: what limits precision? Phase noise ~ 1/SNR per sample, reduced by averaging over a window; and the conversion frequency, which must be the local \(\hat{f}_c\).
I6Derive the attenuation-induced downshift for a Gaussian pulse. Which \(\sigma\), which \(z\), which units? · onsite
Amplitude spectrum \(G(f) = \exp\bigl(-(f-f_0)^2/2\sigma_f^2\bigr)\); round trip to one-way depth \(z\) multiplies by \(\exp(-2\alpha_0 f z)\) with \(\alpha_0\) the amplitude coefficient in Np/cm/MHz. Completing the square: center \(f_0 - 2\alpha_0\sigma_f^2 z\), same \(\sigma_f\). With the power-spectrum \(\sigma_p = \sigma_f/\sqrt{2}\) the same result reads \(f_0 - 4\alpha_0\sigma_p^2 z\). 5 MHz, 60 %, 0.5 dB/cm/MHz (0.0576 Np): −0.187 MHz/cm. Follow-up: turn this into an estimator: \(\hat{\alpha}_0 = -\text{slope}/(2\sigma_f^2)\), with \(\sigma_f\) measured from a reference; then list the confounders (diffraction, non-Gaussian spectra, frequency-dependent backscatter, LPF clipping of the shifted band).
I7For a CUTE-type sound-speed method that measures phase differences between steered IQ images, what biases the delay maps and how do you fix each? · research
(1) Conversion frequency: \(\Delta\tau = -\Delta\varphi/(2\pi f)\), and \(f\) must be the local echo frequency, which decreases with depth (−0.19 MHz/cm in the example above); using \(f_d\) or nominal \(f_0\) scales delays by \(\hat{f}_c/f\). Fix: estimate \(\hat{f}_c\) from fast-time \(R(1)\). (2) Phase wrapping when \(\lvert\Delta\tau\rvert > 1/(2f_c)\): limit the angle step or unwrap. (3) Low-coherence pixels (low envelope, decorrelation between angles) give random phase: weight by \(\lvert R\rvert\) or coherence. (4) Wrong or inconsistent IQ bookkeeping (sign, \(f_d\) in the rotation) between beamforming and phase estimation. Follow-up: why does a constant \(\Delta\) not matter for the phase difference itself? Both images carry the same residual rotation at the same pixel, so it cancels in the product.
I8You are designing a ULM acquisition: 5 MHz probe, 128 channels, 40 mm depth, 10 kHz plane-wave rate. RF or IQ channel data? Quantify. · research
Round trip to 40 mm: 2 × 40/1.54 = 52 µs. RF at 40 MHz, 16-bit: 2078 samples × 128 × 2 B = 532 kB per transmit, 5.3 GB/s at 10 kHz. IQ at 10 MHz: 520 complex = 1040 reals, 266 kB per transmit, 2.7 GB/s, the same as RF at 20 MHz (4 samples per wavelength). IQ is the natural choice: downstream SVD clutter filtering and Doppler need IQ anyway, and IQ DAS with rotation is accurate at 10 MHz (−0.4 dB with linear interpolation). Follow-up: what do you lose? The band outside the LPF; if you want harmonics or full-band spectral analysis, keep RF.
I9How would you estimate the local center frequency robustly from IQ, and why not average the instantaneous frequency? · onsite
Average \(R(1) = \sum \mathrm{IQ}^*[k]\,\mathrm{IQ}[k+1]\) over a depth window and over lines, then \(\hat{f}_c = f_d + (f_{s,\mathrm{IQ}}/2\pi)\arg R(1)\) (08.20). This is equivalent to the \(A^2\)-weighted mean of \(f_i\), which equals the spectral centroid; unweighted \(f_i\) is dominated by spikes at envelope nulls. In the MATLAB lab, 64 lines and 4 mm windows give 0.02 MHz rms scatter and recover \(\alpha_0\) within 3 %. Follow-up: what aliasing limit? \(\lvert f_c - f_d\rvert < f_{s,\mathrm{IQ}}/2\).
I10Your system samples at exactly \(4f_0\) and forms IQ by taking I and Q from neighboring samples. When is that acceptable? · onsite
When the echo is narrowband relative to the \(1/(4f_0)\) offset between the I and Q samples (50 ns at 5 MHz). For a 60 % Gaussian echo the pairwise IQ has 18 % rms complex error averaged over echo positions; averaging the two Q neighbors ([¼ ½ ¼] after the 1, −j, −1, +j mixer) brings it to about 4 %, and a longer interpolating LPF lower still. Follow-up: what does the 3-tap filter do spectrally? \(\cos^2(\pi f/f_s)\): a null at \(f_s/2 = 2f_0\), where the image sits, and 17 % droop at the 2.7 MHz band edge.
I11Why is IQ delay-and-sum exact even though it uses a single frequency \(f_d\) for the rotation? · research
Because the rotation is not a model of the echo; it undoes the demodulator. \(\mathrm{IQ}_n(t) = x_{a,n}(t)\,e^{-j2\pi f_d t}\) exactly, so \(\mathrm{IQ}_n(\tau)\,e^{j2\pi f_d\tau} = x_{a,n}(\tau)\) for any \(\tau\) and all frequencies in the band. The single frequency is the one that was removed, not an approximation of the echo spectrum. Error sources are only interpolation of the baseband and LPF imperfections. Follow-up: what if \(f_d\) varied with depth during demodulation? Rotate by the accumulated demodulator phase \(\theta_d(\tau) = 2\pi\int_0^\tau f_d\,dt\), not \(f_d(\tau)\cdot\tau\).
I12Your phase-based displacement estimates are 15 % low at 4 cm depth in liver and correct near the surface. What is the likely cause? · research
The phase-to-displacement conversion uses the nominal \(f_0\), but the echo is downshifted: at 5 MHz, 60 %, 0.5 dB/cm/MHz, \(f_c(4\ \mathrm{cm}) = 5 - 0.75 = 4.25\) MHz. The phase is \(-2\pi(4.25\ \mathrm{MHz})\tau\), so \(\hat{\tau} = \tau \times 4.25/5 = 0.85\tau\): 15 % low. Fix: convert with the local \(\hat{f}_c\) from fast-time \(R(1)\) (Loupas-style 2-D autocorrelation does this jointly). Follow-up: could the demodulation LPF make it worse? Yes: if the shifted band is clipped, \(\hat{f}_c\) is biased toward \(f_d\).
08.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.
m08_analytic_signal.mFFT-mask analytic signal, envelope, phase and fi of an RF echo; Bedrosian, edge wrap and noisy fi · change f0, bw, t0, noise · runtime ≈ 1 s
%% m08_analytic_signal.m - FFT-based analytic signal of an RF echo
% Purpose : build the analytic signal x_a = x + j*H{x} with an FFT mask
% (no hilbert()), then take envelope |x_a|, instantaneous phase
% and instantaneous frequency of a Gaussian-modulated RF echo.
% Compares the envelope with the true Gaussian envelope, checks the
% mean-frequency identity, and shows three failure modes:
% Bedrosian (low f0, very broadband), circular edge wrap, and
% instantaneous-frequency spikes where the envelope is small.
% Change : f0, bw (fractional -6 dB bandwidth), t0 (echo position), noise.
% Expect : envelope rms error ~0.00 % at 5 MHz / 60 %; ~2.5 % of the lobe
% below 0 Hz at 2 MHz / 120 % (envelope error of a few %);
% A^2-weighted mean f_i = spectral centroid; truncated echo leaks
% energy to the far end of the record.
clear; close all; rng(1);
%% Parameters
fs = 40; % sampling rate (MHz = samples per us)
N = 1024; % record length (samples) -> 25.6 us
t = (0:N-1)'/fs; % time (us)
f0 = 5; % center frequency (MHz)
bw = 0.60; % -6 dB fractional bandwidth of the amplitude spectrum
t0 = 12.8; % echo center (us)
sigf = bw*f0/(2*sqrt(2*log(2))); % std of Gaussian amplitude spectrum (MHz)
sigt = 1/(2*pi*sigf); % std of Gaussian envelope (us)
% FFT analytic-signal mask: DC and Nyquist x1, positive bins x2, negative x0
h = zeros(N,1); h(1) = 1; h(2:N/2) = 2; h(N/2+1) = 1;
analytic = @(x) ifft(fft(x).*h); % works column-wise on an N-row matrix
% Gaussian-modulated echo and its true envelope
echo = @(f0_, sigt_, tc) exp(-(t-tc).^2/(2*sigt_^2)).*cos(2*pi*f0_*(t-tc));
genv = @(sigt_, tc) exp(-(t-tc).^2/(2*sigt_^2));
%% 1. Envelope, phase, instantaneous frequency of a clean echo
x = echo(f0, sigt, t0);
xa = analytic(x);
A = abs(xa); % envelope
phi = unwrap(angle(xa)); % instantaneous phase (rad)
fi = [diff(phi); NaN]*fs/(2*pi); % instantaneous frequency (MHz)
tr = genv(sigt, t0);
m = tr > 0.01;
errEnv = 100*sqrt(sum((A(m)-tr(m)).^2)/sum(tr(m).^2));
fprintf('Re{x_a} - x max abs = %.2e (should be ~1e-16)\n', max(abs(real(xa)-x)));
fprintf('Envelope rms error (5 MHz, 60 %%) = %.4f %%\n', errEnv);
% Mean-frequency identity: A^2-weighted f_i vs centroid of |X_a|^2
w = A(1:end-1).*A(2:end);
fiW = sum(w.*fi(1:end-1))/sum(w);
Xa = fft(xa); fk = (0:N-1)'*fs/N; pos = 1:N/2;
fC = sum(fk(pos).*abs(Xa(pos)).^2)/sum(abs(Xa(pos)).^2);
fprintf('A^2-weighted mean f_i = %.4f MHz\n', fiW);
fprintf('Centroid of |X_a(f)|^2 = %.4f MHz\n', fC);
fprintf('Theory (Gaussian, no clip) = %.4f MHz\n', f0);
%% 2. Bedrosian: low f0, very broadband pulse
f0b = 2; bwb = 1.2; sigfb = bwb*f0b/(2*sqrt(2*log(2))); sigtb = 1/(2*pi*sigfb);
xb = echo(f0b, sigtb, t0); Ab = abs(analytic(xb)); trb = genv(sigtb, t0);
mb = trb > 0.01;
errB = 100*sqrt(sum((Ab(mb)-trb(mb)).^2)/sum(trb(mb).^2));
spill = 0.5*erfc(f0b/sigfb/sqrt(2)); % Phi(-f0/sigma_f)
fprintf('\n2 MHz, 120 %%: sigma_f = %.3f MHz, f0/sigma_f = %.2f\n', sigfb, f0b/sigfb);
fprintf(' lobe fraction below 0 Hz = %.2f %%\n', 100*spill);
fprintf(' envelope rms error = %.2f %%\n', errB);
%% 3. Circular edge artifact: echo cut by the start of the record
tE = 0.15; % echo center near the record start (us)
xe = echo(f0, sigt, tE); Ae = abs(analytic(xe));
xe0 = [xe; zeros(N,1)]; % same echo, zero-padded to 2N
h2 = zeros(2*N,1); h2(1) = 1; h2(2:N) = 2; h2(N+1) = 1;
Ae0 = abs(ifft(fft(xe0).*h2)); Ae0 = Ae0(1:N);
tail = t > t(end) - 0.5; % last 0.5 us of the record
fprintf('\nTruncated echo at %.2f us: max |x_a| in last 0.5 us = %.3f (no pad)\n', tE, max(Ae(tail)));
fprintf(' with 2x zero pad = %.3f\n', max(Ae0(tail)));
%% 4. Noise: instantaneous frequency is erratic where the envelope is small
xn = x + 0.05*randn(N,1);
xan = analytic(xn); An = abs(xan);
fin = angle(xan(2:end).*conj(xan(1:end-1)))*fs/(2*pi); % no unwrap needed
wn = An(1:end-1).*An(2:end);
win = abs(t(1:end-1) - t0) < 3*sigt; % echo window, +-3 sigma_t
strong = An(1:end-1) > 0.5*max(An);
weak = An(1:end-1) < 0.1*max(An);
fprintf('\nNoise 5 %% of peak:\n');
fprintf(' std f_i where A > 50 %% max = %.2f MHz\n', std(fin(strong)));
fprintf(' std f_i where A < 10 %% max = %.2f MHz\n', std(fin(weak)));
fprintf(' A^2-weighted mean f_i, echo window = %.3f MHz (clean: %.3f)\n', sum(wn(win).*fin(win))/sum(wn(win)), fiW);
fprintf(' unweighted mean f_i, echo window = %.3f MHz\n', mean(fin(win)));
%% Figures
figure('Position', [100 100 900 650]);
zw = abs(t - t0) < 6*sigt;
subplot(2,2,1);
plot(t(zw), x(zw), 'b', t(zw), imag(xa(zw)), 'r', t(zw), A(zw), 'k', t(zw), -A(zw), 'k', ...
t(zw), tr(zw), 'g--');
xlabel('time (\mus)'); ylabel('amplitude');
title(sprintf('x (b), Hilbert (r), |x_a| (k), Gaussian (g--); err %.3f %%', errEnv));
subplot(2,2,2);
fz = zw & A > 0.1*max(A);
plot(t(fz), fi(fz), 'b.-'); hold on; plot(t(fz), f0+0*t(fz), 'k--');
xlabel('time (\mus)'); ylabel('f_i (MHz)'); title('instantaneous frequency (A > 10 % max)');
ylim([f0-1 f0+1]);
subplot(2,2,3);
plot(t, 20*log10(Ae+1e-12), 'r', t, 20*log10(Ae0+1e-12), 'k--'); ylim([-80 5]);
xlabel('time (\mus)'); ylabel('|x_a| (dB)');
title('truncated echo: no pad (r) vs 2x zero pad (k--)'); legend('no pad', 'zero pad');
subplot(2,2,4);
plot(t(zw), fin(min(find(zw), N-1)), 'b.'); hold on;
plot(t(zw), An(zw)/max(An)*4 + f0 - 2, 'k');
xlabel('time (\mus)'); ylabel('f_i (MHz)'); ylim([-5 15]);
title('noisy echo: f_i (dots) and scaled envelope (k)');
What to look for: \(\mathrm{Re}\{x_a\} - x \approx 2\times10^{-16}\); envelope error 0.004 % at 5 MHz, 60 %; \(A^2\)-weighted \(f_i\) and the \(\lvert X_a\rvert^2\) centroid both 5.000 MHz; at 2 MHz, 120 %: 2.49 % of the lobe below 0 Hz and 2.5 % envelope error; truncated echo leaves 0.095 at the far end without padding, 0.000 with; with 5 % noise, std \(f_i\) 0.54 MHz where \(A\) > 50 % and 9.1 MHz where \(A\) < 10 %.
m08_iq_demodulation.mMix, hand-designed windowed-sinc LPF, decimate; spectra per stage, envelope error, aliasing vs M, fd mismatch · change fd, fc, L, M · runtime ≈ 2 s
%% m08_iq_demodulation.m - IQ demodulation: mix, low-pass, decimate
% Purpose : build the IQ chain by hand: complex mixing with exp(-j2*pi*fd*t),
% a windowed-sinc FIR low-pass you design yourself (Hamming window
% written out, no toolbox), and decimation by M. Shows the spectrum
% at each stage, compares 2|IQ| with the analytic-signal envelope,
% sweeps M to show where aliasing starts, and shows what a
% mixing-frequency mismatch does (residual rotation, band clipping).
% Change : f0, bw, fd (demodulation frequency), fc (LPF cutoff), L (taps), M.
% Expect : M = 4 (10 MHz IQ): envelope error ~0.6 % (set by the LPF);
% reconstruction error stays at that floor to M = 5 (8 MHz) and
% grows from M = 6; ~4 % at M = 8 (5 MHz), ~12 % at M = 9-10.
% fd = f0 + 1 MHz rotates IQ at -1 MHz with the envelope intact
% (1.7 %) while fc = 4 MHz passes the shifted band; fc = 2 MHz
% clips it (envelope error ~24 %, rotation estimate biased).
clear; close all; rng(1);
%% Parameters
fs = 40; N = 2048; t = (0:N-1)'/fs; % MHz, samples, us
f0 = 5; bw = 0.6; % echo center (MHz), -6 dB frac. BW
sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf);
fd = 5; % demodulation frequency (MHz)
fc = 4; % LPF cutoff (MHz)
L = 49; % FIR length (odd -> integer delay)
M = 4; % decimation factor
% RF line: a few echoes of different phase plus weak noise
tc = [8 11 11.4 14 20]; amp = [1 0.6 -0.8 0.4 0.7];
x = zeros(N,1);
for k = 1:numel(tc)
x = x + amp(k)*exp(-(t-tc(k)).^2/(2*sigt^2)).*cos(2*pi*f0*(t-tc(k)) + k);
end
x = x + 0.002*randn(N,1);
% Reference: analytic signal by FFT mask
h = zeros(N,1); h(1) = 1; h(2:N/2) = 2; h(N/2+1) = 1;
xa = ifft(fft(x).*h);
% Windowed-sinc low-pass (Hamming window written out), unit DC gain
% sinc without the toolbox: sin(pi x)/(pi x), with sinc(0) = 1
sinc_ = @(u) (sin(pi*u) + (u==0)) ./ (pi*u + (u==0));
designLPF = @(fc_, L_) ( (2*fc_/fs)*sinc_((-(L_-1)/2:(L_-1)/2)'*2*fc_/fs) ...
.* (0.54 - 0.46*cos(2*pi*(0:L_-1)'/(L_-1))) );
normDC = @(b) b/sum(b);
lpf = @(v, b) conv(v, b, 'same'); % zero-phase alignment for odd L
% Full chain as one anonymous function: returns IQ at fs (before decimation)
iqchain = @(x_, fd_, fc_) 2*lpf(x_.*exp(-1j*2*pi*fd_*t), normDC(designLPF(fc_, L)));
%% 1. The chain at nominal settings
b = normDC(designLPF(fc, L));
xm = x.*exp(-1j*2*pi*fd*t); % stage 2: mixed
xl = 2*lpf(xm, b); % stage 3: filtered (x2 restores amplitude)
iq = xl(1:M:end); fsIQ = fs/M; % stage 4: decimated
ref = xa.*exp(-1j*2*pi*fd*t); % ideal IQ = analytic signal, carrier removed
% Envelope error at the kept samples
idx = (1:M:N)'; sel = abs(xa(idx)) > 0.05;
eKept = 100*sqrt(sum((abs(iq(sel)) - abs(xa(idx(sel)))).^2)/sum(abs(xa(idx(sel))).^2));
% Band-limited reconstruction back to 40 MHz (what an interpolator sees)
Nd = numel(iq); rec = interpft(iq, Nd*M);
eRec = 100*sqrt(sum(abs(abs(rec(1:N)) - abs(xa)).^2.*(abs(xa)>0.05))/sum(abs(xa).^2.*(abs(xa)>0.05)));
% Complex error vs ideal IQ (includes phase)
ws = abs(xa) > 0.05; % evaluate where the echoes are (noise floor excluded)
eCplx = 100*norm(xl(ws) - ref(ws))/norm(ref(ws));
fprintf('fd = %.1f MHz, fc = %.1f MHz, L = %d, M = %d (IQ rate %.1f MHz)\n', fd, fc, L, M, fsIQ);
fprintf(' 2|IQ| vs |x_a| at kept samples = %.3f %% rms\n', eKept);
fprintf(' band-limited reconstruction error = %.3f %% rms\n', eRec);
fprintf(' complex IQ error before decimation = %.3f %% rms\n', eCplx);
fprintf(' -20 dB half-width of the echo band = %.2f MHz (needs fs/M > %.1f MHz)\n', ...
sigf*sqrt(2*log(10)), 2*sigf*sqrt(2*log(10)));
%% 2. Sweep the decimation factor
Ms = 1:12; eM = zeros(size(Ms));
for i = 1:numel(Ms)
q = xl(1:Ms(i):end); nd = numel(q);
r = interpft(q, nd*Ms(i)); r = r(1:min(N, nd*Ms(i)));
rr = xa(1:numel(r)); w = abs(rr) > 0.05;
eM(i) = 100*sqrt(sum((abs(r(w)) - abs(rr(w))).^2)/sum(abs(rr(w)).^2));
end
fprintf('\nM : '); fprintf('%6d', Ms); fprintf('\nfsIQ: '); fprintf('%6.1f', fs./Ms);
fprintf('\nerr (%%): '); fprintf('%6.2f', eM); fprintf('\n');
%% 3. Mixing-frequency mismatch
fdm = 6; % demodulate 1 MHz above f0
for fcm = [4 2]
xlm = iqchain(x, fdm, fcm);
ephm = 100*sqrt(sum((abs(xlm) - abs(xa)).^2.*(abs(xa)>0.05))/sum(abs(xa).^2.*(abs(xa)>0.05)));
% residual rotation rate from the lag-one product inside the first echo
w8 = abs(t - tc(1)) < 2*sigt;
rot = angle(sum(conj(xlm(find(w8)-1)).*xlm(w8)))*fs/(2*pi);
fprintf('fd = %.1f MHz, fc = %.1f MHz: envelope error %.2f %%, residual rotation %.3f MHz (expect %.1f)\n', ...
fdm, fcm, ephm, rot, f0 - fdm);
end
%% Figures: spectra at each stage, envelopes in time
F = @(v) 20*log10(abs(fftshift(fft(v)))/max(abs(fft(x))) + 1e-12);
fax = ((0:N-1)' - N/2)*fs/N;
Hb = 20*log10(abs(fftshift(fft(b, N))) + 1e-12);
figure('Position', [80 80 900 700]);
subplot(3,2,1); plot(fax, F(x)); axis([-20 20 -80 5]); title('1 RF (real): lobes at \pm f_0'); ylabel('dB');
subplot(3,2,2); plot(fax, F(xm)); axis([-20 20 -80 5]); title('2 mixed: baseband + image at -(f_0+f_d)');
subplot(3,2,3); plot(fax, F(xl/2), 'b', fax, Hb, 'k--'); axis([-20 20 -80 5]);
title(sprintf('3 after LPF (f_c = %g MHz, %d taps), filter dashed', fc, L)); ylabel('dB');
Nq = numel(iq); faq = ((0:Nq-1)' - Nq/2)*fsIQ/Nq;
subplot(3,2,4); plot(faq, 20*log10(abs(fftshift(fft(iq)))/max(abs(fft(iq))) + 1e-12));
axis([-fsIQ/2 fsIQ/2 -80 5]); title(sprintf('4 decimated by %d: span \\pm%.2f MHz', M, fsIQ/2)); xlabel('frequency (MHz)');
subplot(3,2,5); tw = t > 6 & t < 16;
plot(t(tw), abs(xa(tw)), 'k', t(idx(idx>=find(tw,1) & idx<=find(tw,1,'last'))), ...
abs(iq(idx>=find(tw,1) & idx<=find(tw,1,'last'))), 'ro');
xlabel('time (\mus)'); ylabel('envelope'); title('|x_a| (k) and 2|IQ| at kept samples (o)');
subplot(3,2,6); plot(Ms, eM, 'bo-'); xlabel('decimation M'); ylabel('reconstruction error (%)');
title('band-limited reconstruction error vs M'); grid on;
figure('Position', [100 100 800 300]);
xlm = iqchain(x, fdm, 4);
plot(t(tw), real(xl(tw)), 'b', t(tw), real(xlm(tw)), 'r', t(tw), abs(xlm(tw)), 'k');
xlabel('time (\mus)'); ylabel('amplitude');
title('I with f_d = f_0 (b) and f_d = f_0 + 1 MHz (r): residual 1 MHz rotation, same envelope (k)');
What to look for: at \(M\) = 4, \(2\lvert\mathrm{IQ}\rvert\) matches \(\lvert x_a\rvert\) to 0.56 % at the kept samples and 0.60 % after band-limited reconstruction. The \(M\) sweep stays at 0.6 % to \(M\) = 5, then 1.4, 3.8, 3.9, 12 % for \(M\) = 6 to 9. With \(f_d\) = 6 MHz the residual rotation measures −0.99 MHz and the envelope error is 1.7 % (\(f_c\) = 4 MHz) versus 24 % with \(f_c\) = 2 MHz.
m08_iq_beamforming.m32-element plane-wave point target: RF DAS vs IQ DAS with and without rotation, lateral profiles, loss vs depth · change Nel, z0, M, fd · runtime ≈ 20 s
%% m08_iq_beamforming.m - DAS on RF vs DAS on IQ, with and without phase rotation
% Purpose : simulate plane-wave (0 deg) channel data for one point target on
% a linear array, demodulate each channel to IQ (mix, FIR low-pass,
% decimate by 4), and beamform three ways:
% (a) RF DAS (RF upsampled x4, linear interpolation) + analytic
% signal along depth -> reference envelope
% (b) IQ DAS: interpolate IQ at tau_n, rotate by exp(+j2*pi*fd*tau_n)
% (c) IQ DAS without the rotation
% Prints peak loss of (b) and (c) re (a), compares (c) with the
% narrowband sum |mean(exp(-j2*pi*fd*tau_n))| and the Fresnel
% integral, and sweeps depth.
% Change : Nel (32-64), z0, f0, M (IQ decimation), fd (try fd ~= f0).
% Expect : 32 el, 5 MHz, z = 20 mm: (b) -0.4 dB re (a) (LPF droop plus
% linear interpolation of IQ at 10 MHz; -0.1 dB with M = 2);
% (c) -12.4 dB image peak, -13.3 dB on the target (narrowband
% prediction -13.4 dB, Fresnel -12.5 dB); the loss oscillates with
% depth and is still -6.5 dB at 60 mm. Runtime ~20 s in Octave.
clear; close all; rng(1);
%% Parameters
c = 1.54; % speed of sound (mm/us)
f0 = 5; bw = 0.6; % MHz, -6 dB fractional bandwidth
fd = f0; % demodulation frequency (MHz)
fs = 40; M = 4; fsIQ = fs/M; % RF and IQ sampling rates (MHz)
Nel = 32; pitch = 0.3; % elements, pitch (mm)
z0 = 20; % target depth (mm)
sigf = bw*f0/(2*sqrt(2*log(2))); sigt = 1/(2*pi*sigf); % MHz, us
xel = ((0:Nel-1) - (Nel-1)/2)*pitch; % element x (mm), row
% FIR low-pass for the IQ chain (49-tap Hamming windowed sinc, fc = 4 MHz)
L = 49; fcut = 4; k = (-(L-1)/2:(L-1)/2)';
u = 2*fcut/fs*k; sc = (sin(pi*u) + (u==0))./(pi*u + (u==0));
b = (2*fcut/fs)*sc.*(0.54 - 0.46*cos(2*pi*(0:L-1)'/(L-1))); b = b/sum(b);
% Channel data (columns = elements). Plane-wave transmit: t = (z + r)/c
simChannels = @(zt, xt, Nt) exp(-((0:Nt-1)'/fs - (zt + sqrt(zt^2 + (xel - xt).^2))/c).^2/(2*sigt^2)) ...
.* cos(2*pi*f0*((0:Nt-1)'/fs - (zt + sqrt(zt^2 + (xel - xt).^2))/c));
%% Beamform: point target at (0, z0)
Nt = ceil((2*z0 + 12)/c*fs) + 64; Nt = Nt + mod(Nt, 2*M);
t = (0:Nt-1)'/fs;
rf = simChannels(z0, 0, Nt);
iqf = 2*conv2(rf.*exp(-1j*2*pi*fd*t), b, 'same'); % mix + LPF (all channels)
iq = iqf(1:M:end, :); tIQ = (0:size(iq,1)-1)'/fsIQ; % decimate
rfU = interpft(rf, 4*Nt); tU = (0:4*Nt-1)'/(4*fs); % RF upsampled to 160 MHz
xg = linspace(-3, 3, 121); % lateral grid (mm)
zg = z0 + linspace(-1, 1, 128)'; % depth grid (mm), 15.7 um step
A = zeros(numel(zg), numel(xg)); B = A; C = A;
for ix = 1:numel(xg)
tau = (zg + sqrt(zg.^2 + (xg(ix) - xel).^2))/c; % Nz x Nel delays (us)
sa = zeros(numel(zg),1); sb = sa; scc = sa;
for n = 1:Nel
sa = sa + interp1(tU, rfU(:,n), tau(:,n), 'linear', 0);
v = interp1(tIQ, real(iq(:,n)), tau(:,n), 'linear', 0) + ...
1j*interp1(tIQ, imag(iq(:,n)), tau(:,n), 'linear', 0);
sb = sb + v.*exp(1j*2*pi*fd*tau(:,n));
scc = scc + v;
end
A(:,ix) = sa; B(:,ix) = sb; C(:,ix) = scc;
end
% analytic signal of the RF-beamformed column along depth (zero-padded FFT mask)
Np = 2*numel(zg); h = zeros(Np,1); h(1) = 1; h(2:Np/2) = 2; h(Np/2+1) = 1;
Aa = ifft(fft([A; zeros(Np-numel(zg), numel(xg))]).*repmat(h, 1, numel(xg)));
Aenv = abs(Aa(1:numel(zg), :));
pkA = max(Aenv(:)); pkB = max(abs(B(:))); pkC = max(abs(C(:)));
taun = (z0 + sqrt(z0^2 + xel.^2))/c;
pred = 20*log10(abs(mean(exp(-1j*2*pi*fd*taun))));
lam = c/f0; psi = pi*max(xel)^2/(lam*z0);
s = linspace(0, 1, 4001); fr = 20*log10(abs(trapz(s, exp(-1j*psi*s.^2))));
fprintf('Nel = %d, pitch %.2f mm, f0 = fd = %.1f MHz, z0 = %d mm, IQ rate %.1f MHz\n', Nel, pitch, f0, z0, fsIQ);
fprintf(' (b) IQ + rotation peak re (a) = %+.2f dB\n', 20*log10(pkB/pkA));
fprintf(' (c) IQ, no rotation peak re (a) = %+.2f dB\n', 20*log10(pkC/pkA));
fprintf(' narrowband prediction |mean(exp(-j2pi fd tau_n))| = %+.2f dB\n', pred);
fprintf(' Fresnel psi = %.2f rad, |int_0^1 exp(-j psi s^2) ds| = %+.2f dB\n', psi, fr);
% lateral -6 dB width of (a) and (b)
latA = max(Aenv)/pkA; latB = max(abs(B))/pkA; latC = max(abs(C))/pkA;
w6 = @(p) (xg(2)-xg(1))*sum(p >= 0.5*max(p));
fprintf(' lateral -6 dB width: (a) %.2f mm, (b) %.2f mm, (c) %.2f mm (lambda F# = %.2f mm)\n', ...
w6(latA), w6(latB), w6(latC), lam*z0/(Nel*pitch));
%% Depth sweep of the no-rotation loss (peak at the target only)
zs = 10:5:60; lossC = zeros(size(zs)); lossB = lossC; predS = lossC; frS = lossC;
for i = 1:numel(zs)
zt = zs(i); Nt2 = ceil((2*zt + 12)/c*fs) + 64; Nt2 = Nt2 + mod(Nt2, 2*M); t2 = (0:Nt2-1)'/fs;
r2 = simChannels(zt, 0, Nt2);
q2 = 2*conv2(r2.*exp(-1j*2*pi*fd*t2), b, 'same'); q2 = q2(1:M:end,:); tq = (0:size(q2,1)-1)'/fsIQ;
zz = zt + linspace(-0.3, 0.3, 61)'; tau = (zz + sqrt(zz.^2 + xel.^2))/c;
sb = zeros(size(zz)); sc2 = sb; sr = sb;
for n = 1:Nel
v = interp1(tq, real(q2(:,n)), tau(:,n), 'linear', 0) + 1j*interp1(tq, imag(q2(:,n)), tau(:,n), 'linear', 0);
sb = sb + v.*exp(1j*2*pi*fd*tau(:,n)); sc2 = sc2 + v;
end
tn = (zt + sqrt(zt^2 + xel.^2))/c;
lossC(i) = 20*log10(max(abs(sc2))/max(abs(sb)));
predS(i) = 20*log10(abs(mean(exp(-1j*2*pi*fd*tn))));
ps = pi*max(xel)^2/(lam*zt); frS(i) = 20*log10(abs(trapz(s, exp(-1j*ps*s.^2))));
end
fprintf('\nDepth sweep (on-axis peak of (c) re (b)):\n z (mm) '); fprintf('%7d', zs);
fprintf('\n (c)-(b) '); fprintf('%7.1f', lossC);
fprintf('\n NB sum '); fprintf('%7.1f', predS);
fprintf('\n Fresnel '); fprintf('%7.1f', frS); fprintf('\n');
%% Figures
figure('Position', [80 80 950 600]);
dbI = @(I) max(20*log10(abs(I)/pkA + 1e-12), -40);
subplot(2,3,1); imagesc(xg, zg, dbI(Aenv)); caxis([-40 0]); colormap(gray);
xlabel('x (mm)'); ylabel('z (mm)'); title('(a) RF DAS');
subplot(2,3,2); imagesc(xg, zg, dbI(B)); caxis([-40 0]); xlabel('x (mm)'); title('(b) IQ + rotation');
subplot(2,3,3); imagesc(xg, zg, dbI(C)); caxis([-40 0]); xlabel('x (mm)'); title('(c) IQ, no rotation'); colorbar;
subplot(2,3,4:5);
plot(xg, 20*log10(latA), 'k', xg, 20*log10(latB), 'b--', xg, 20*log10(latC), 'r'); ylim([-50 3]);
xlabel('lateral x (mm)'); ylabel('dB re (a) peak'); legend('(a) RF', '(b) IQ + rot', '(c) IQ no rot');
title('lateral profile (max over depth)'); grid on;
subplot(2,3,6); plot(zs, lossC, 'ro-', zs, predS, 'k.--', zs, frS, 'b:'); grid on;
xlabel('depth (mm)'); ylabel('loss (dB)'); legend('sim', 'NB sum', 'Fresnel', 'Location', 'southeast');
title('no-rotation peak loss vs depth');
What to look for: (b) −0.42 dB re (a) (−0.14 dB with \(M\) = 2); (c) −12.4 dB at the image peak and −13.3 dB on target, against −13.45 dB narrowband and −12.48 dB Fresnel; lateral −6 dB width 0.75 mm for (a) and (b), smeared to 6 mm for (c). The depth sweep tracks the narrowband sum within 0.2 dB from 10 to 60 mm (−14.1 dB at 10 mm, −6.5 dB at 60 mm).
m08_attenuation_downshift.mGaussian spectrum × e−2α0fz: centroid vs −2α0σf²z, then a toy spectral-shift attenuation estimate from speckle IQ · change alpha_dB, f0, bw, Nl · runtime ≈ 1 s
%% m08_attenuation_downshift.m - attenuation downshift and a toy QUS estimate
% Purpose : (1) multiply a Gaussian amplitude spectrum (center f0, std sigma_f)
% by the round-trip attenuation exp(-2*alpha0*f*z) and compare the
% measured power-spectrum centroid with f0 - 2*alpha0*sigma_f^2*z.
% alpha0 is the AMPLITUDE attenuation slope in Np/cm/MHz, z is the
% ONE-WAY depth in cm, sigma_f is the std of the AMPLITUDE spectrum
% (= sqrt(2) x the std of the power spectrum).
% (2) toy QUS: simulate speckle RF lines whose pulse downshifts with
% depth, demodulate to IQ, estimate the local center frequency with
% the fast-time lag-one autocorrelation, fit the slope versus depth
% and invert it for alpha0.
% Change : alpha_dB (dB/cm/MHz), f0, bw, window length, number of lines.
% Expect : slope -0.187 MHz/cm at 5 MHz / 60 % / 0.5 dB/cm/MHz;
% part (1) matches theory to < 0.001 MHz; part (2) recovers
% alpha0 = 0.49 dB/cm/MHz (true 0.50; a few % spread from speckle
% and from the LPF/decimation choices).
clear; close all; rng(1);
%% Parameters
f0 = 5; bw = 0.6; % MHz, -6 dB fractional BW (amplitude)
sigf = bw*f0/(2*sqrt(2*log(2))); % std of the amplitude spectrum (MHz)
alpha_dB = 0.5; % dB/cm/MHz (liver-like)
alpha = alpha_dB/(20*log10(exp(1))); % Np/cm/MHz (1 Np = 8.686 dB)
fprintf('sigma_f = %.3f MHz (power-spectrum sigma_p = %.3f MHz), alpha0 = %.4f Np/cm/MHz\n', ...
sigf, sigf/sqrt(2), alpha);
%% 1. Spectral model: centroid versus depth
f = linspace(0, 20, 4001)'; % MHz (positive frequencies only)
z = 0:0.5:8; % one-way depth (cm)
G = exp(-(f - f0).^2/(2*sigf^2)); % transmitted amplitude spectrum
S = G .* exp(-2*alpha*f*z); % received amplitude spectra (columns)
P = S.^2; % power spectra
fcP = sum(f.*P)./sum(P); % power-spectrum centroid (MHz)
fcA = sum(f.*S)./sum(S); % amplitude-spectrum centroid, for comparison
fcT = f0 - 2*alpha*sigf^2*z; % theory
sdP = sqrt(sum((f - fcP).^2.*P)./sum(P)); % power-spectrum std (should stay sigma_f/sqrt 2)
fprintf('\nTheory slope -2*alpha0*sigma_f^2 = %.4f MHz/cm\n', -2*alpha*sigf^2);
fprintf(' z (cm) fc power fc ampl theory sigma_p\n');
for i = 1:2:numel(z)
fprintf('%6.1f %7.3f %7.3f %7.3f %7.3f\n', z(i), fcP(i), fcA(i), fcT(i), sdP(i));
end
fprintf('max |fc power - theory| = %.2e MHz\n', max(abs(fcP - fcT)));
lossdB = 20*log10(exp(-2*alpha*f0*z(end)));
fprintf('Round-trip loss at f0, z = %.0f cm: %.1f dB (= 2 x %.1f dB/cm/MHz x %.0f MHz x %.0f cm)\n', ...
z(end), lossdB, alpha_dB, f0, z(end));
%% 2. Toy QUS: speckle lines with a depth-dependent pulse
c = 1.54; fs = 40; M = 4; fsIQ = fs/M; fd = f0; % mm/us, MHz
zmax = 80; Nt = round(2*zmax/c*fs); t = (0:Nt-1)'/fs; zt = c*t/2; % depth axis (mm)
Nl = 64; % RF lines (independent speckle)
blk = 2; % pulse is held constant over 2 mm blocks
sigt = 1/(2*pi*sigf); tp = (-round(5*sigt*fs):round(5*sigt*fs))'/fs;
nb = ceil(zmax/blk);
RF = zeros(Nt, Nl);
for ib = 1:nb
rows = find(zt >= (ib-1)*blk & zt < ib*blk);
zc = ((ib-0.5)*blk)/10; % block center (cm)
fcz = f0 - 2*alpha*sigf^2*zc; % local center frequency (closed form)
pulse = exp(-tp.^2/(2*sigt^2)).*cos(2*pi*fcz*tp); % TGC: amplitude loss removed
scat = zeros(Nt, Nl); scat(rows, :) = randn(numel(rows), Nl);
RF = RF + conv2(scat, pulse, 'same');
end
% IQ chain: mix, 49-tap Hamming windowed-sinc LPF (fc = 4.5 MHz), decimate by M
L = 49; fcut = 4.5; k = (-(L-1)/2:(L-1)/2)'; u = 2*fcut/fs*k;
b = (2*fcut/fs)*(sin(pi*u) + (u==0))./(pi*u + (u==0)).*(0.54 - 0.46*cos(2*pi*(0:L-1)'/(L-1)));
b = b/sum(b);
IQ = 2*conv2(RF.*exp(-1j*2*pi*fd*t), b, 'same'); IQ = IQ(1:M:end, :); zIQ = zt(1:M:end);
% Lag-one autocorrelation in 4 mm windows, averaged over lines
win = 4; zc = (win/2 + 5):2:(zmax - win/2 - 5); fhat = zeros(size(zc));
for i = 1:numel(zc)
r = find(abs(zIQ - zc(i)) < win/2);
R1 = sum(sum(conj(IQ(r(1:end-1), :)).*IQ(r(2:end), :)));
fhat(i) = fd + fsIQ/(2*pi)*angle(R1);
end
p = polyfit(zc/10, fhat, 1); % slope in MHz/cm
alphaHat = -p(1)/(2*sigf^2);
fprintf('\nToy QUS (%d lines, %g mm windows, IQ at %.0f MHz):\n', Nl, win, fsIQ);
fprintf(' fitted slope = %.4f MHz/cm (theory %.4f)\n', p(1), -2*alpha*sigf^2);
fprintf(' alpha0 estimate = %.4f Np/cm/MHz = %.3f dB/cm/MHz (true %.3f)\n', ...
alphaHat, alphaHat*20*log10(exp(1)), alpha_dB);
fprintf(' rms scatter of f_hat about the fit = %.3f MHz\n', std(fhat - polyval(p, zc/10)));
%% Figures
figure('Position', [80 80 950 360]);
subplot(1,3,1); plot(f, S(:, 1:4:end)./max(S(:,1))); xlim([0 12]);
xlabel('frequency (MHz)'); ylabel('amplitude (norm. to z = 0)');
title('received spectra, z = 0, 2, 4, 6, 8 cm');
subplot(1,3,2); plot(z, fcP, 'bo', z, fcT, 'k-', z, fcA, 'r.'); grid on;
xlabel('one-way depth z (cm)'); ylabel('center frequency (MHz)');
legend('power centroid', 'f_0 - 2\alpha_0\sigma_f^2 z', 'amplitude centroid', 'Location', 'southwest');
title('downshift: model check');
subplot(1,3,3); plot(zc/10, fhat, 'b.', zc/10, polyval(p, zc/10), 'r-', zc/10, f0 - 2*alpha*sigf^2*zc/10, 'k--');
grid on; xlabel('depth (cm)'); ylabel('f_c estimate (MHz)');
title(sprintf('toy QUS: \\alpha_0 = %.2f dB/cm/MHz', alphaHat*20*log10(exp(1))));
legend('R(1) estimate', 'linear fit', 'theory', 'Location', 'southwest');
What to look for: slope −0.1869 MHz/cm; power centroid matches theory within \(2\times10^{-4}\) MHz to 8 cm while the power-spectrum \(\sigma_p\) stays 0.901 MHz; round-trip loss 40 dB at 8 cm; the toy QUS fit gives −0.182 MHz/cm and \(\hat{\alpha}_0\) = 0.487 dB/cm/MHz (true 0.5).
08.13Self-check
Q1A 5 MHz pulse has 80 % (−6 dB) fractional bandwidth and a Gaussian spectrum. What is the minimum complex IQ sampling rate, and how does it compare with real RF sampling?
\(B_{-6} = 0.8 \times 5 = 4\) MHz, so the bare minimum (−6 dB support) is \(f_{s,\mathrm{IQ}} > 4\) MHz with \(f_d = f_0\). A Gaussian is not band-limited: \(\sigma_f = 4/2.355 = 1.70\) MHz, the −20 dB full width is \(2\sigma_f\sqrt{2\ln 10}\) = 7.3 MHz and the −40 dB width is 10.3 MHz. A practical choice is 8 to 10 MHz complex (16 to 20 million real numbers per second). Real sampling needs \(f_s > 2f_{\max}\) with \(f_{\max} = 5 + 3.65 = 8.65\) MHz at −20 dB, i.e. ≥ 17.3 MHz, and 30 to 40 MHz in practice for good interpolation. The IQ rate is set by bandwidth, the RF rate by the highest frequency.
Q2An echo is delayed by 10 ns at 5 MHz. What phase rotation does its IQ signal show, and what axial displacement does that represent?
\(\Delta\varphi = -2\pi f_0\tau\) \(= -2\pi \times 5\times10^{6} \times 10\times10^{-9}\) \(= -0.314\) rad = −18°. Displacement \(d = c\tau/2 = 1540 \times 10\times10^{-9}/2\) = 7.7 µm. The envelope shift of 10 ns is 0.4 samples at 40 MHz and invisible in the magnitude; the phase measures it directly.
Q3Your system demodulates with \(f_d\) = 5.3 MHz but the echo is centered at 5.0 MHz. What happens to the baseband spectrum, the envelope and phase-based estimates?
\(\Delta = f_0 - f_d = -0.3\) MHz, so the baseband lobe is centered at −0.3 MHz instead of 0 and the image at \(-(f_0+f_d) = -10.3\) MHz. In time, \(\mathrm{IQ}(t) = b(t)\,e^{-j2\pi\cdot 0.3 t}\): a residual rotation of 0.3 cycles per µs (one turn every 3.3 µs, ≈ 2.6 mm of depth). The envelope is unchanged as long as the LPF passes \(\lvert\Delta\rvert + B/2\) and \(f_{s,\mathrm{IQ}} > B + 0.6\) MHz. Slow-time Doppler phase is unaffected (the residual rotation is identical in every frame at a given depth), but converting phase to delay or velocity with \(f_d\) instead of the true 5.0 MHz biases it: the measured \(\Delta\varphi = -2\pi(5.0\ \mathrm{MHz})\tau\), so \(\hat{\tau} = -\Delta\varphi/(2\pi f_d) = \tau \times 5.0/5.3\), an underestimate of 5.7 %. An IQ beamformer must still rotate with \(f_d\) = 5.3 MHz, the frequency actually removed.
Q4Why does the FFT Hilbert mask keep bins 0 and \(N/2\) at weight 1 instead of doubling or zeroing them?
For real \(x\), \(X[0]\) and \(X[N/2]\) are real and each is its own mirror (\(-0 = 0\) and \(-N/2 \equiv N/2 \bmod N\)). They are shared between the positive and negative halves, so they are neither redundant nor to be doubled. Weight 1 preserves \(\mathrm{Re}\{x_a\} = x\) exactly, which is the defining property of the analytic signal. Doubling DC would add a constant to \(x\); zeroing it would remove the mean from the real part. Bandpass RF has negligible energy there, so the effect is small in practice, but the mask should be right.
Q5You compute the envelope of a 2 MHz, 120 %-bandwidth pulse with hilbert() and it does not match the Gaussian you transmitted. Is the code wrong?
No. \(\sigma_f = 2.4/2.355 = 1.02\) MHz and \(f_0/\sigma_f = 1.96\), so about \(\Phi(-1.96) \approx 2.5\) % of the Gaussian lobe lies below 0 Hz. The analytic signal discards it and reflects the conjugate of the negative lobe's tail onto the positive side, so \(x_a \ne a(t)\,e^{j2\pi f_0 t}\) (Bedrosian's condition fails). The "true envelope" of such a pulse is not uniquely defined; \(\lvert x_a\rvert\) is the standard, well-defined choice. Demo 08.2 shows the error directly.
Q6An IQ beamformer forgets the phase rotation. Estimate the peak loss at \(z\) = 20 mm for 32 elements, 0.3 mm pitch, 5 MHz.
\(x_{\max} = 15.5 \times 0.3 = 4.65\) mm, \(\lambda\) = 0.308 mm. \(\psi = \pi x_{\max}^2/(\lambda z)\) \(= \pi \times 21.6/(0.308 \times 20)\) = 11.0 rad, so the aperture spans ≈ 1.75 cycles of quadratic phase (3.5 Fresnel zones). \(\bigl\lvert\int_0^1 e^{-j\psi s^2}ds\bigr\rvert \approx 0.24\) → −12.5 dB; the exact 32-element sum gives −13.4 dB (Demo 08.5 reproduces this). The PSF also broadens because the uncompensated phase acts like a strong defocus.
Q7At 5 MHz, 60 % bandwidth, the fast-time centroid falls by 0.20 MHz per cm of depth. What attenuation coefficient does this imply?
\(\sigma_f = 3/2.355 = 1.274\) MHz, \(\sigma_f^2 = 1.623\) MHz². \(\alpha_0 = -\text{slope}/(2\sigma_f^2)\) = 0.20/3.246 = 0.0616 Np/cm/MHz = 0.535 dB/cm/MHz (× 8.686). This is in the range of normal-to-fatty liver. Using the power-spectrum \(\sigma_p\) = 0.901 MHz with the form \(4\alpha_0\sigma_p^2\) gives the same answer; mixing the two \(\sigma\) definitions is the classic factor-of-2 error.
Q8You decimate IQ by \(M\) = 8 from 40 MHz and compare \(2\lvert\mathrm{IQ}[k]\rvert\) with \(\lvert x_a[8k]\rvert\): they match to 0.1 %. Is the decimation safe?
Not necessarily. Decimation picks samples without altering them, so the kept samples always match the reference regardless of aliasing. Aliasing shows up only in what you compute between samples. The test is whether the LPF output is band-limited to \(\pm f_{s,\mathrm{IQ}}/2\) = ±2.5 MHz: for a 5 MHz, 60 % pulse the −20 dB half-width is 2.7 MHz, so it is not, and band-limited reconstruction (and therefore beamforming interpolation) is wrong by several percent. Demo 08.3 reports both errors.
08.14Figures from the literature
Measured data to hold this module's models against. Images are not reproduced; open the linked figure. Figure numbers were checked against the linked full text.
Centroid downshift with depth, measured in human liver
K. J. Parker, R. M. Lerner, R. C. Waag, "Comparison of techniques for in vivo attenuation measurements," IEEE Trans. Biomed. Eng. 35(12):1064–1068, 1988. Full text (author site). Fig. 5(a), 5(c), 5(d) and Fig. 2.
Fig. 5 is representative liver data: (a) averaged spectra from four depths, (c) center frequency versus depth from the centroid and zero-crossing estimators with linear fits, (d) spectral standard deviation versus depth. Fig. 2 compares the centroid-shift attenuation with amplitude decay across 49 subjects. Look at 5(c) as the measured version of (08.19): a linear downshift whose slope, divided by \(2\sigma_f^2\), is the attenuation. Then look at 5(d): the Gaussian model says \(\sigma_f\) stays constant with depth, so any trend there is the departure from that model that a real transducer spectrum and diffraction produce. The paper reports centroid-shift values of 0.18 to 1.55 dB/MHz·cm with about 0.2 dB/MHz·cm intraliver scatter, wider than amplitude decay (0.25 to 0.90), the noise price of a spectral estimator.
Echo phase between steered complex images, in vivo sound speed
P. Stähli, M. Kuriakose, M. Frenz, M. Jaeger, "Forward model for quantitative pulse-echo speed-of-sound imaging," arXiv:1902.10639v2 (2019); published as "Improved forward model for quantitative pulse-echo speed-of-sound imaging," Ultrasonics 108, 2020. Fig. 2(c), Figs. 4–6 and Fig. 7 (arXiv numbering).
The method is this module's phase-to-delay step applied to images: each RF signal is made analytic with the Hilbert transform, complex images are beamformed per transmit angle, and the local echo phase shift is the angle of the kernel-averaged Hermitian product of two images (2 × 2 mm kernel, 2° angle steps to avoid phase aliasing). Fig. 2(c) sketches how the measured phase shift follows the shift of the echo-envelope centroid; Figs. 4–6 are phantom sound-speed maps for three layer geometries and three forward models; Fig. 7 is a volunteer's abdominal wall and liver. Read the angle-step choice as the \(\lvert\tau\rvert < 1/(2f_0)\) wrap limit of (08.13), and the kernel as the averaging that tames phase noise in low-envelope pixels.
IQ delay-and-sum feeding vector Doppler on a rotating disk
V. Perrot, M. Polichetti, F. Varray, D. Garcia, "So you think you can DAS? A viewpoint on delay-and-sum beamforming," Ultrasonics 111:106309, 2021. doi:10.1016/j.ultras.2020.106309 (preprint arXiv:2007.11960). Fig. 12 and Fig. 14 (arXiv v3 numbering).
Fig. 12 is experimental: vector Doppler in a rotating disk from 17 unsteered plane waves, with the I/Q signals DAS-beamformed at ±15° receive angles, and the median velocity error plotted against f-number. Fig. 14 shows how a sparse DAS matrix produces interpolated I/Q along the receive hyperbolas. This is (08.15) in production: interpolation of the slow baseband plus the phase rotation, followed by slow-time phase processing that only works if every pixel's phase bookkeeping is consistent.
Re-beamforming the same channel data across sound speeds, in vivo
V. Perrot, M. Polichetti, F. Varray, D. Garcia, Ultrasonics 111:106309, 2021 (as above). doi:10.1016/j.ultras.2020.106309. Fig. 8, Fig. 10 and Fig. 11 (arXiv v3 numbering).
Fig. 8 uses the PICMUS CIRS phantom data to show contrast and lateral resolution versus the assumed speed of sound and f-number; Figs. 10 and 11 apply the same sweep to eight carotids and report the speed of sound that maximizes an image-quality metric. Each point of such a sweep is a full re-beamforming of one channel data set: the case Section 08.7 gives for IQ beamforming, where the rotation \(e^{j2\pi f_d\tau_n(c)}\) makes every re-delay exact at a low sample rate and where a wrong \(c\) shows up as phase incoherence across the aperture.
08.15Go deeper
- A. V. Oppenheim, R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson, 2010. Ch. 12 (Discrete Hilbert transforms): analytic sequences, the discrete Hilbert kernel, FIR Hilbert transformers.
- R. G. Lyons, Understanding Digital Signal Processing, 3rd ed., Prentice Hall, 2011. Ch. 8 (Quadrature signals) and Ch. 9 (The discrete Hilbert transform): the clearest treatment of complex mixing and negative frequency.
- S. L. Marple, "Computing the discrete-time 'analytic' signal via FFT," IEEE Trans. Signal Processing 47(9):2600–2603, 1999. The DC/Nyquist mask.
- E. Bedrosian, "A product theorem for Hilbert transforms," Proc. IEEE 51(5):868–869, 1963.
- T. L. Szabo, Diagnostic Ultrasound Imaging: Inside Out, 2nd ed., Academic Press, 2014. Attenuation and its effect on pulse spectra; signal processing chain of a scanner.
- C. Kasai, K. Namekawa, A. Koyano, R. Omoto, "Real-time two-dimensional blood flow imaging using an autocorrelation technique," IEEE Trans. Sonics Ultrason. SU-32(3):458–464, 1985.
- T. Loupas, J. T. Powers, R. W. Gill, "An axial velocity estimator for ultrasound blood flow imaging, based on a full evaluation of the Doppler equation by means of a two-dimensional autocorrelation approach," IEEE Trans. UFFC 42(4):672–688, 1995.
- R. Kuc, "Estimating acoustic attenuation from reflected ultrasound signals: comparison of spectral-shift and spectral-difference approaches," IEEE Trans. Acoust., Speech, Signal Process. 32(1):1–6, 1984.
- V. Perrot, M. Polichetti, F. Varray, D. Garcia, "So you think you can DAS? A viewpoint on delay-and-sum beamforming," Ultrasonics 111:106309, 2021. Includes IQ beamforming and the phase-rotation term.
- M. Jaeger, G. Held, S. Peeters, S. Preisser, M. Grünig, M. Frenz, "Computed ultrasound tomography in echo mode for imaging speed of sound using pulse-echo sonography: proof of principle," Ultrasound Med. Biol. 41(1):235–250, 2015.