Doppler and flow imaging
Every flow mode measures one thing: how the phase of a backscattered echo changes from one transmit to the next. Spectral, colour, power, ultrafast and vector Doppler differ in how they sample that slow-time signal, how they remove the tissue that sits 40 dB above it, and which estimator they apply. The same slow-time signal, filtered by SVD, is what ULM tracks and what sound-speed and aberration estimators must not corrupt.
You will be able to
- Derive the Doppler equation, the pulsed-Doppler Nyquist limit and the range–velocity product \(v_{\max} z_{\max} = c^{2}/(8 f_0\cos\theta)\), and quote the numbers at 3, 5 and 15 MHz
- Predict the shape of a Doppler spectrum from the velocity profile, the gate, the elevation beam and the aperture (geometric broadening)
- Write the Kasai and Loupas autocorrelators, their bias and variance, and the noise bias of the variance (turbulence) estimate
- Design a wall filter (regression or SVD) for a given ensemble and tissue motion, and explain why ultrafast power Doppler gains sensitivity as \(\sqrt{N}\)
- Estimate 2D velocity vectors from multi-angle Doppler or transverse oscillation, and predict the lateral precision penalty \(\cot\gamma\)
Assumes
- Module 06 · Eigendecomposition & SVD: SVD clutter filtering; Module 08 · Analytic signal & IQ demodulation: IQ demodulation
- Module 09 · Digital filter design: regression wall filters; Module 11 · Correlation & time-delay estimation: the lag-one phase estimator
- Module 23 · Beam patterns, k-space and compounding: plane-wave coherent compounding; Module 28 · Microbubble physics, contrast imaging and the ULM signal source: microbubbles as flow tracers
29.1The Doppler equation and why the angle dominates the error budget
A scatterer moving with velocity component \(v\cos\theta\) along the beam shortens (or lengthens) the round-trip path by \(2v\cos\theta\) per second. The echo phase \(\varphi = -4\pi f_0 z/c\) therefore rotates at the Doppler frequency
- \(f_D\)
- Doppler shift, Hz (positive for flow toward the probe)
- \(v\)
- scatterer speed, m/s
- \(f_0\)
- transmit (centre) frequency, Hz
- \(\theta\)
- angle between the beam and the velocity vector
- \(c\)
- speed of sound, 1540 m/s
The equation is exact to first order in \(v/c\); the second-order term \((v/c)^{2}\) is of order \(10^{-6}\) and never matters. What matters is \(\cos\theta\). The scanner measures \(f_D\), the operator sets \(\theta\), and the velocity is \(v = c f_D/(2 f_0\cos\theta)\). An angle error \(\delta\theta\) produces
- \(\theta\)
- beam-to-flow angle, rad
- \(\delta v,\ \delta\theta\)
- velocity error (m/s) and angle error (rad)
- \(v_{\mathrm{est}},\ v\)
- estimated and true speed, m/s
- \(\theta_{\mathrm{true}},\ \theta_{\mathrm{set}}\)
- true beam-to-flow angle and the angle-correction setting on the scanner, rad
CW versus PW. A continuous-wave probe transmits on one half of the aperture and receives on the other. The Doppler signal comes from everything in the overlap of the two beams: no range resolution, but no sampling either, so no aliasing (the limit is the analog bandwidth). Pulsed-wave Doppler transmits a short burst every \(T = 1/\mathrm{PRF}\) and samples the demodulated echo once per pulse at the delay \(2z/c\): this range gate buys depth selectivity at the price of sampling the Doppler signal at the PRF. Everything in 29.2 follows from that sampling.
29.2Pulsed Doppler: the PRF, the Nyquist limit and the range–velocity trade-off
Between two pulses the scatterer moves \(vT\cos\theta\) toward the probe and the gated IQ sample rotates by \(\Delta\varphi = 4\pi f_0 v T\cos\theta/c = 2\pi f_D T\). The slow-time sequence \(z_n\) is a complex sinusoid at \(f_D\) sampled at the PRF, so it is unambiguous only while \(\lvert\Delta\varphi\rvert < \pi\):
- \(\mathrm{PRF}\)
- pulse repetition frequency (slow-time sampling rate), Hz
- \(v_{\mathrm{Nyq}}\)
- Nyquist (aliasing) velocity, m/s
- \(f_D,\ v,\ f_0,\ \theta,\ c\)
- as in (29.1)
The next pulse cannot be fired until the echo from the gate has returned, \(T \ge 2z/c\). Combining \(\mathrm{PRF} \le c/(2z_{\max})\) with (29.3) removes the PRF:
- \(z_{\max}\)
- deepest unambiguous gate depth, m
- \(v_{\max}\)
- largest unaliased speed at that depth, m/s
- \(c,\ f_0,\ \theta\)
- as in (29.1)
- \(\mathrm{PRF}\)
- pulse repetition frequency, Hz
When the flow exceeds \(v_{\mathrm{Nyq}}\) the spectrum wraps: the peak of a systolic waveform reappears at the opposite edge of the display. Three remedies, in order of preference: shift the baseline (the band is still one PRF wide but can be placed asymmetrically, up to nearly \(2v_{\mathrm{Nyq}}\) in one direction), raise the PRF (until \(z_{\max}\) reaches the gate; beyond that "high-PRF" mode accepts several gates in range), or lower \(f_0\) or increase \(\theta\) (the last one trades aliasing for angle error).
Spectral Doppler: aliasing, baseline shift and the depth limit
Approx. modelParabolic flow insonated over the whole lumen (flat spectrum from 0 to \(f_{D,\max}(t)\)), exponential speckle per time–frequency cell, illustrative pulsatile waveform with a 1.0 m/s systolic peak at 60 beats per minute. The wrapping of the spectrum by the PRF is exact.29.3Spectral Doppler: what shapes the spectrum
The Doppler power spectrum is, to first order, the histogram of Doppler frequencies of the scatterers in the sample volume, weighted by how strongly the beam and gate insonate each one:
- \(W(\mathbf{r})\)
- two-way beam and gate weight at position \(\mathbf{r}\) in the vessel cross-section
- \(v(\mathbf{r})\)
- local speed along the vessel axis, m/s
- \(S(f)\)
- Doppler power spectrum at frequency \(f\), arbitrary units
- \(\propto\)
- proportional to
- \(\delta\)
- Dirac delta function
- \(f_0,\ \theta,\ c\)
- as in (29.1)
- \(d^{2}\mathbf{r}\)
- area element of the vessel cross-section, m²
Three things change that shape. Profile: plug flow gives a narrow line at \(f_{D,\max}\); a blunted profile, \(1 - (r/a)^{4}\), puts the mean at \(\tfrac23 f_{D,\max}\). Gate and elevation beam: a small gate on the axis keeps only the fast core; a thin elevation beam samples one diameter rather than the disk, which tilts the spectrum toward \(f_{D,\max}\) (mean \(\tfrac23\) for the full diameter). The aperture: each element sees the scatterer at a slightly different angle, so a single velocity produces a band of frequencies (geometric or "intrinsic" spectral broadening). With transmit and receive angles spread over \(\theta \pm \alpha\), where \(\alpha = \arctan(1/(2F_{\#}))\),
- \(\alpha\)
- aperture half-angle seen from the gate, rad
- \(F_{\#}\)
- F-number (focal depth / aperture width)
- \(f,\ f_D\)
- frequency of a spectral component and the Doppler shift of (29.1), Hz
- \(\theta\)
- beam-to-flow angle, rad
The transit-time picture gives the same answer: a scatterer crossing a beam of width \(W\) at speed \(v\sin\theta\) is seen for \(W/(v\sin\theta)\) seconds, and a signal of that duration has a bandwidth of order \(v\sin\theta/W\). Newhouse and colleagues showed in 1980 that geometric and transit-time broadening are one phenomenon described two ways. The point for a velocity estimator is the same: broadening scales with \(v\), so it is a fixed fraction of \(f_D\) set by the aperture and \(\tan\theta\).
The spectrogram is a short-time Fourier transform. The display is \(\lvert\mathrm{STFT}\rvert^{2}\) of the slow-time IQ sequence with a window of \(N_w\) samples:
- \(\Delta t\)
- observation time of one spectrum, s
- \(N_w\)
- number of pulses in the analysis window
- \(\mathrm{PRF}\)
- pulse repetition frequency, Hz
- \(\Delta f_{-6\,\mathrm{dB}}\)
- spectral resolution (−6 dB width) with a Hann window, Hz
The Doppler spectrum of a vessel: profile, gate, elevation beam and aperture
SimulationMonte Carlo histogram of 60 000 scatterers placed uniformly in the gated part of the lumen, each with transmit and receive angles drawn uniformly in \(\theta \pm \alpha\); same model as m29_doppler_spectrum_profile.m. The geometry panel is drawn to scale for a 6 mm vessel.29.4Colour flow: short ensembles and the autocorrelation estimators
Colour flow imaging cannot afford a spectrum per pixel. It fires a short ensemble (packet) of \(N\) = 6–16 pulses along each line, and estimates three numbers per pixel from the slow-time sequence: the mean Doppler frequency, the power and the spectral width. The estimator, due to Kasai and colleagues (1985), is the phase of the lag-one autocorrelation, the same estimator you met for displacement in Module 11:
- \(z_n\)
- complex (IQ) sample of one pixel at pulse \(n\), after the wall filter
- \(N\)
- ensemble length (pulses per pixel)
- \(\hat R(1)\)
- lag-one slow-time autocorrelation estimate
- \(\arg\)
- four-quadrant phase, in \((-\pi, \pi]\)
- \(c,\ f_0\)
- as in (29.1)
- \(\mathrm{PRF}\)
- pulse repetition frequency, Hz
- \(\hat v\)
- estimated axial velocity, m/s
- \((\cdot)^{*}\)
- complex conjugate
- \(n\)
- pulse (slow-time) index
The hidden frequency. (29.8) uses \(f_0\), the demodulation frequency. The received echo is centred lower, because attenuation downshifts the spectrum (Module 27, frequency-dependent attenuation) and each scatterer weights frequencies differently. Since the Doppler phase per pulse is \(4\pi f v T/c\) at every RF frequency \(f\), the measured phase reflects the mean RF frequency \(\bar f_{\mathrm{RF}}\): using \(f_0\) biases \(\hat v\) by \(\bar f_{\mathrm{RF}}/f_0\). Loupas, Powers and Gill (1995) estimate \(\bar f_{\mathrm{RF}}\) from the lag-one correlation along depth on the same data block:
- \(z_{m,n}\)
- IQ sample at depth index \(m\) and pulse \(n\)
- \(f_s\)
- IQ sampling rate along depth, Hz
- \(f_{\mathrm{dem}}\)
- demodulation frequency, Hz
- \(\hat R(0,1),\ \hat R(1,0)\)
- lag one in slow time, lag one in depth
- \(c,\ \mathrm{PRF}\)
- as in (29.1), (29.3)
- \(\hat v\)
- estimated velocity, m/s
- \(m\prime,\ n\prime\)
- lag in depth and slow time, samples
- \((\cdot)^{*}\)
- complex conjugate
Spectral width and turbulence. The magnitude of the normalized correlation measures how fast the signal decorrelates. For a Gaussian spectrum of standard deviation \(\sigma_f\), \(\lvert R(T)\rvert/R(0) = e^{-(2\pi\sigma_f T)^{2}/2}\), so
- \(\hat\sigma_\omega\)
- spectral width (angular), rad/s; \(\sigma_f = \sigma_\omega/2\pi\)
- \(T\)
- pulse interval 1/PRF, s
- \(\mathrm{SNR}\)
- signal-to-noise power ratio per sample
- \(\hat R(0),\ \hat R(1)\)
- lag-zero (power) and lag-one slow-time autocorrelation, as in (29.8)
- \(\lvert\cdot\rvert\)
- magnitude
The Kasai estimator on a short ensemble: N, noise, spectral width and aliasing
Simulation400 independent ensembles per setting: complex Gaussian slow-time process with a Gaussian Doppler spectrum (white noise shaped by FFT, aliased into one PRF) plus white noise; Kasai and variance estimators applied to the first \(N\) samples. \(f_0\) = 5 MHz, PRF 4 kHz, \(v_{\mathrm{Nyq}}\) = 0.308 m/s, \(\theta\) = 0.29.5Wall filters: removing tissue 40 dB above the blood
Before any estimator, the tissue echo (the vessel wall, myocardium, liver) has to go. It is 20–40 dB stronger than blood in most organs and more near strong walls, and it moves slowly: breathing, pulsation and probe motion give tissue Doppler of a few to tens of hertz. Left in place it drags \(\arg\hat R(1)\) toward the tissue velocity and floods the power map. The wall filter is a high-pass filter in slow time, but with \(N\) = 8–16 samples there is no room for a long FIR or for IIR transients. The standard answer is the polynomial regression filter (Module 09): project out the best-fitting polynomial of degree \(p\) over the ensemble,
- \(\mathbf{z},\ \mathbf{y}\)
- ensemble before and after filtering, \(N\times 1\) complex
- \(\mathbf{Q}\)
- \(N\times(p+1)\) orthonormal basis of polynomials of degree \(\le p\) in slow time
- \(p\)
- regression order (0 = mean removal)
- \(\mathbf{I}\)
- identity matrix
- \((\cdot)^{T}\)
- transpose
- \(f_{-3\,\mathrm{dB}}\)
- −3 dB cut-off of the wall filter, Hz
- \(\mathrm{PRF},\ N\)
- pulse repetition frequency (Hz) and ensemble length
The filter has \(N - p - 1\) degrees of freedom left, so it also removes \((p+1)/N\) of the white noise, and every degree of \(p\) raises the minimum detectable velocity. The design choice is a trade between residual clutter (too low a cutoff: flash artefacts and velocities pulled toward zero) and lost slow flow (too high a cutoff: holes near the vessel walls where the parabolic profile is slow). The demo shows this on a simulated colour-flow frame. Sections 29.6 and the SVD filter remove the trade-off by using space as well as time.
Colour-flow map with a regression wall filter
SimulationSame model as m29_color_flow_wallfilter.m, smaller grid: inclined vessel (radius 3 mm, \(v_{\max}\) = 0.25 m/s, \(\theta\) = 60°, parabolic), tissue speckle 40 dB above blood moving axially at \(v_t\), blood-to-noise 15 dB; each signal a sum of 8 random phasors around its mean Doppler (blood σ = 40 Hz, tissue σ = 3 Hz). \(f_0\) = 5 MHz, PRF 2 kHz, \(N\) = 12. Detection: \(\hat R(0)\) more than 3 dB above the expected filtered noise power.29.6Power Doppler, ultrafast Doppler and functional ultrasound
Power Doppler displays \(\hat R(0)\), the energy left after the wall filter, and ignores the phase. It has no aliasing and no angle dependence of its value (only of whether the flow survives the wall filter), and the power is proportional to the number of moving red cells in the resolution cell: a blood-volume map rather than a velocity map. Its sensitivity is set by how well the power estimate separates a weak vessel from the noise floor, and that improves with the number of independent samples.
Ultrafast Doppler (Bercoff et al., IEEE TUFFC 2011) replaces line-by-line packets with tilted plane waves that insonify the whole field, coherently compounded into frames at 0.5–2 kHz (Figure 29.B). Every pixel then has a continuous slow-time record: the ensemble is set by the acquisition time, not by the number of lines, and \(N\) = 200–400 is routine. Two things follow. First, the wall filter can use space: tissue motion is spatially coherent (the whole organ moves together), blood is not, so the Casorati matrix of pixels × time separates them by rank (Demené et al., IEEE TMI 2015; Module 06):
- \(\mathbf{S}\)
- Casorati matrix, \(P\) pixels × \(N\) frames, complex IQ
- \(\sigma_k,\ \mathbf{u}_k,\ \mathbf{v}_k\)
- singular values, spatial and temporal singular vectors
- \(K\)
- clutter rank: number of components removed
- \(\mathbf{V}_K\)
- first \(K\) temporal singular vectors, \(N\times K\)
- \(\mathbf{S}_f\)
- filtered (clutter-removed) Casorati matrix
- \(\mathbf{I}\)
- identity matrix
- \((\cdot)^{H}\)
- conjugate transpose
- \(N\)
- number of frames (ensemble length)
Second, the power estimate gets \(N - K\) degrees of freedom. For complex Gaussian noise of power \(P_n\) the background power Doppler value has mean \(P_n (N-K)/N\) and standard deviation \(P_n\sqrt{N-K}/N\); a vessel adds \(P_b (N-K)/N\). The contrast-to-noise ratio of the power map is therefore
- \(\hat P\)
- power Doppler value of a pixel
- \(P_b,\ P_n\)
- blood and noise power per sample (after compounding)
- \(N - K\)
- independent noise samples left after the filter
- \(\mathbf{r},\ n\)
- pixel position and frame index
- \(s_f(\mathbf{r}, n)\)
- IQ sample after clutter filtering
- \(N,\ K\)
- number of frames and clutter rank removed, as in (29.12)
- \(\bar P_{\mathrm{vessel}},\ \bar P_{\mathrm{bg}}\)
- mean power in the vessel and background regions
- \(\mathrm{std}\)
- standard deviation over the background region
- \(\mathrm{CNR}\)
- contrast-to-noise ratio, dimensionless
The assumption hidden in \(\sqrt{N-K}\) is that blood decorrelates within the ensemble. Blood in a 1 mm/s microvessel at 15 MHz has a Doppler of 19.5 Hz and stays correlated for tens of milliseconds; over a 25-frame ensemble it looks as static as tissue and part of it is projected out with the clutter. The script shows this: at \(N\) = 25 the 1 mm/s vessel reaches a CNR of 0.41 against 1.2 predicted, at \(N\) = 400 it is within 16 % of theory.
Functional ultrasound (Macé et al., Nature Methods 2011) uses exactly this sensitivity: ultrafast plane-wave power Doppler images transient changes of cerebral blood volume, including in vessels too small for conventional Doppler, and those changes during a stimulus (neurovascular coupling) are the functional signal. The first paper mapped whisker-evoked activity in the rat cortex and thalamus and followed the spread of epileptic seizures. ULM (Module 31) takes the same SVD-filtered data and, instead of integrating power, localizes the individual microbubbles in it.
Ultrafast power Doppler: SVD clutter filter and sensitivity versus ensemble length
SimulationSame model as m29_ultrafast_power_doppler.m on a 32 × 32 pixel grid (100 µm pixels): tissue speckle 40 dB above blood with non-rigid axial motion (50 µm at 3.3 Hz, ×1.5 at the bottom); four 200 µm vessels at 1, 3, 10 and 30 mm/s (Gaussian Doppler spectra); noise 6 dB above blood; 15 MHz, 1 kHz frames. Top singular components by subspace iteration on \(\mathbf{S}^{H}\mathbf{S}\) with Rayleigh–Ritz, JS.29.7Vector flow: measuring the component the beam cannot see
Every estimator so far measures \(v\cos\theta\). Flow perpendicular to the beam gives zero Doppler, and the operator's angle correction assumes the flow direction is known (parallel to the vessel wall), which fails in bifurcations, stenoses, aneurysms and the heart. Vector flow imaging measures two (or three) components. Four families, reviewed by Jensen, Nikolov, Yu and Garcia (IEEE TUFFC 2016):
- Multi-angle (cross-beam) Doppler. Transmit and/or receive from \(K \ge 2\) directions; each pair \((\alpha_k, \beta_k)\) measures a different projection of the velocity:
\[ f_k = -\frac{f_0}{c}\Big[(\sin\alpha_k + \sin\beta_k)\,v_x + (\cos\alpha_k + \cos\beta_k)\,v_z\Big], \qquad \hat{\mathbf{v}} = (\mathbf{A}^{T}\mathbf{A})^{-1}\mathbf{A}^{T}\hat{\mathbf{f}}, \quad \mathrm{Cov}(\hat{\mathbf{v}}) = \sigma_f^{2}(\mathbf{A}^{T}\mathbf{A})^{-1} \](29.14)
- \(\alpha_k,\ \beta_k\)
- transmit and receive steering angles of pair \(k\) (from the array normal)
- \(v_x,\ v_z\)
- lateral and axial velocity (z downward, so \(v_z < 0\) is toward the probe), m/s
- \(\mathbf{A}\)
- \(K\times 2\) matrix of the bracketed coefficients times \(-f_0/c\)
- \(\sigma_f\)
- standard deviation of each Doppler estimate, Hz
- \(f_0,\ c\)
- as in (29.1)
- \(f_k\)
- Doppler shift measured for pair \(k\), Hz
- \(k\)
- index of the transmit/receive angle pair, \(1\dots K\)
- \(\hat{\mathbf{v}},\ \hat{\mathbf{f}}\)
- estimated velocity vector \((v_x, v_z)\) and the vector of measured shifts \(f_k\)
- \(\mathrm{Cov}\)
- covariance matrix
- \((\cdot)^{T},\ (\cdot)^{-1}\)
- transpose, matrix inverse
- Transverse oscillation (Jensen and Munk 1998; Anderson 1998). A receive apodization with two peaks separated by \(d\) produces a field that oscillates laterally with wavelength
- \(\lambda_x\)
- lateral oscillation wavelength, m
- \(z\)
- depth, m
- \(d\)
- separation of the two apodization peaks, m
- \(\lambda\)
- wavelength, m
- Speckle tracking. 2D block matching between frames (Module 11), angle independent but limited by speckle decorrelation (out-of-plane motion, shear, flow gradients) and by the frame rate; plane-wave imaging made it practical for blood.
- Directional beamforming and model-based methods. Beamform along the flow direction (directional cross-correlation), or combine Doppler with the continuity equation and wall boundary conditions (vector flow mapping in echocardiography).
The common lesson is geometric: the lateral component is measured through an effective aperture angle (\(\gamma\) or \(d/z\)) much smaller than the 180° of the axial "aperture", so its precision is several times worse and degrades with depth. The same \(\cot\gamma\) penalty governs lateral displacement estimation in elastography and the lateral resolution of sound-speed tomography from limited angles (Module 32).
Multi-angle vector Doppler: least squares, covariance and the cot γ penalty
Simulation400 trials: per angle, \(N\) slow-time samples of a complex exponential at \(f_k\) with random phase plus white noise, Kasai estimate, then least-squares vector; ellipses are the theoretical 1σ and 2σ contours \(\sigma_f^{2}(\mathbf{A}^{T}\mathbf{A})^{-1}\) with the measured per-angle \(\sigma_f\). Same model as m29_vector_flow_multiangle.m: 5 MHz, transmit rate 10 kHz shared by the \(K\) angles.29.8Aliasing and unwrapping, turbulence, and the Doppler SNR budget
Recognizing aliasing in a colour map. An aliased velocity jumps from \(+v_{\mathrm{Nyq}}\) to \(-v_{\mathrm{Nyq}}\) between neighbouring pixels: the colour goes from the brightest red straight to the brightest blue. A true flow reversal passes through zero, so it goes through black. The estimate is the true velocity wrapped into one PRF,
- \(m\)
- unknown wrap count per pixel
- \(T_1,\ T_2\)
- two interleaved pulse intervals (staggered PRF), s
- \(v_{\mathrm{ext}}\)
- extended unambiguous velocity of the dual-PRF difference, m/s
- \(v,\ \hat v\)
- true velocity and aliased (wrapped) estimate, m/s
- \(v_{\mathrm{Nyq}}\)
- as in (29.3), m/s
- \(\mathbb{Z}\)
- set of integers
- \(\lambda\)
- wavelength at \(f_0\), m
- \(\lvert\cdot\rvert\)
- magnitude
Spatial unwrapping works when the velocity field is smooth over the pixel grid and the region is connected; it fails at the vessel wall (noise) and in jets narrower than a few pixels. Robust implementations unwrap along paths of high power and low variance, which is the same quality-guided strategy used for phase unwrapping in MRI and interferometry.
Turbulence and spectral broadening. Disturbed flow downstream of a stenosis has a wide range of velocities in the sample volume, so its spectrum widens and the colour "variance" rises. Two other effects widen it too and must be subtracted before calling it turbulence: the geometric broadening of 29.3 (a fixed fraction \(\approx 2\alpha\tan\theta\) of \(f_D\)) and noise (29.10). A velocity gradient across the sample volume (the near-wall shear) also widens the spectrum without any turbulence.
Doppler SNR budget. Blood backscatter is 20–40 dB below the surrounding tissue (Jensen, Estimation of Blood Velocities Using Ultrasound, 1996), so flow imaging lives near the noise floor. What buys SNR, in dB:
| Lever | Gain | Cost |
|---|---|---|
| Ensemble length \(N\) (power map CNR) | \(5\log_{10}N\) in CNR (×2 per 4×) | time per frame; blood must decorrelate |
| Coherent compounding of \(n_a\) plane waves | up to \(10\log_{10}n_a\) in SNR (4.8 dB for 3) | PRF divided by \(n_a\): \(v_{\mathrm{Nyq}}\) drops |
| Longer transmit burst (\(M\) cycles) | energy ∝ \(M\): +3 dB per doubling | axial resolution, gate length |
| Spatial averaging of \(\hat R(1)\) over \(n_s\) pixels | variance of \(\hat v\) ∝ \(1/n_s\) | resolution of the colour map |
| Lower \(f_0\) | less attenuation (≈ 0.5 dB/cm/MHz one way in soft tissue) | lower \(f_D\), wider PSF; blood scattering ∝ \(f^{4}\) falls |
29.9Visual summary
Twelve questions an expert answers in five seconds. Each figure is computed from the formula or model named in its badge.
29.10How experts think about it
A flow expert sees a slow-time signal per pixel, sampled at the PRF, made of a huge, slow, spatially coherent tissue component, a weak, decorrelating blood component and white noise. Every design question becomes: how many samples, at what rate, which subspace to remove, and which moment of the remaining spectrum to estimate. The angle, the Nyquist limit and the noise floor are the three things that go wrong.
The round-trip path changes at \(2v\cos\theta\), so the phase per pulse is \(4\pi f_0 vT\cos\theta/c\) (29.1). At 90° the Doppler shift is zero whatever the speed: a vessel parallel to the array is invisible to axial Doppler unless the beams are steered. 1 m/s at 60° and 5 MHz gives 3.25 kHz.
The correction divides by \(\cos\theta\); its derivative gives \(\tan\theta\,\delta\theta\) (29.2). The exact error \(\cos\theta/\cos(\theta+\delta\theta) - 1\) for a 5° misalignment is 5.7 % at 30°, 18.3 % at 60° and 32 % at 70°. The vessel axis, not the jet direction, is what the operator aligns, so in a stenosis the true angle is unknown anyway.
\(\mathrm{PRF} \le c/(2z)\) and \(v_{\mathrm{Nyq}} = c\,\mathrm{PRF}/(4f_0\cos\theta)\) multiply to (29.4). At 10 cm and 3 MHz the limit is 0.99 m/s; a 4 m/s aortic jet at 10 cm would need \(f_0\) = 0.74 MHz. Baseline shift buys at most a factor 2; beyond that, CW or high-PRF.
The estimator returns a phase in \((-\pi, \pi]\), so a velocity of 1.3 \(v_{\mathrm{Nyq}}\) reads −0.7 \(v_{\mathrm{Nyq}}\). In a colour map the alias boundary joins the two brightest colours; a real reversal passes through black. A baseline shift moves the window (figure: shifted by +0.4 PRF), it does not widen it.
With \(v = v_{\max}(1 - r^{2}/a^{2})\), the area element \(2\pi r\,dr\) is proportional to \(dv\), so every velocity is equally represented (29.5). A thin elevation beam samples one diameter (\(dv \propto u\,du\)) and gives density \(\propto 1/\sqrt{1 - v/v_{\max}}\), rising to the peak; a blunted \(1 - (r/a)^{4}\) profile under a wide beam gives the same density. Volume flow = mean velocity × area only when the whole lumen is uniformly insonated.
Each transmit–receive element pair sees the scatterer at a slightly different angle (29.6). With \(\alpha = \arctan(1/2F_{\#})\) = 14° and \(\theta\) = 60° the edge is at \(\cos 46°/\cos 60°\) = 1.39 \(f_D\). The spectral-envelope "peak velocity" overestimates by up to that factor; vendors calibrate it out, which is why peak velocities differ between scanners.
\(R(T)\) is the Fourier transform of the Doppler spectrum evaluated at one lag. For a Gaussian spectrum its phase is \(2\pi\bar f T\) and its magnitude falls as a Gaussian in \(\sigma_f T\) (29.10). The variance estimate is useful only while \(\sigma_f T\) is small but not negligible: at PRF 4 kHz a 150 Hz width moves \(\lvert R(1)\rvert/R(0)\) by just 3 %, so a small bias in \(R(0)\) (noise) dominates.
The Doppler phase per pulse is \(4\pi f vT/c\) at every frequency in the echo band, so the measured phase is set by the mean RF frequency, which attenuation pulls below \(f_0\) in proportion to depth, attenuation slope and the square of the bandwidth (Module 27). Loupas divides by \(f_{\mathrm{dem}} + (f_s/2\pi)\arg R(1,0)\) (29.9). Figure: m29_kasai_loupas.m, 400 trials each.
Noise adds power to \(R(0)\) and nothing to \(R(1)\) (it is white), so \(\lvert R(1)\rvert/R(0)\) shrinks by \(\mathrm{SNR}/(1+\mathrm{SNR})\) and the width estimate grows (29.10). Figure: apparent versus true width at PRF 4 kHz. At 10 dB SNR the floor alone (true width zero) is 271 Hz; any width below that is unmeasurable. Subtracting an estimate of the noise power from \(R(0)\) removes most of the bias.
Projecting out \(p+1\) polynomials removes a band around DC whose width scales as PRF/N (29.11). Each order adds about 0.41 PRF/N of cutoff and about 20–25 dB of stop-band depth at low frequency (figure, computed). To lower \(v_{\min}\) without losing rejection, lengthen \(N\), not lower \(p\).
An organ moving as a whole gives every pixel the same temporal signature (up to a slowly varying gain), so it occupies a few temporal singular vectors (29.12); blood speckle in different vessels decorrelates independently and spreads over all \(N\). Figure: the Demo 29.5 model at \(N\) = 128, four components 20–55 dB above a flat floor.
The background power estimate is an average of \(N - K\) independent exponential variables: its standard deviation falls as \(1/\sqrt{N-K}\) while the vessel's excess power stays (29.13). Figure: m29_ultrafast_power_doppler.m for the 10 mm/s vessel (dots) against theory (dashed): 1.23, 1.54, 2.61, 3.44, 4.68 for \(N\) = 25 … 400.
A signal at \(f_D\) completes one phase turn in \(\mathrm{PRF}/f_D\) samples. Over a shorter ensemble it is nearly constant and lies in the same temporal subspace as static tissue; SVD projects part of it out and a high-pass removes it entirely. Figure: \(\mathrm{PRF}/f_D\) versus speed at 15 MHz and 1 kHz. The script shows the consequence: the 1 mm/s vessel reaches CNR 0.41 at \(N\) = 25 against 1.2 predicted.
Two symmetric pairs at \(\pm\gamma\) give rows \(-(2f_0/c)[\pm\sin\gamma, \cos\gamma]\); \(\mathbf{A}^{T}\mathbf{A}\) is diagonal with \(2\sin^{2}\gamma\) and \(2\cos^{2}\gamma\) (29.14). With receive fixed at 0° the rows are \([\pm\sin\gamma, 1 + \cos\gamma]\) and the ratio is \(\cot(\gamma/2)\). m29_vector_flow_multiangle.m reproduces 5.66 at 10° (theory 5.67).
The receive field in the focal plane is the Fourier transform of the apodization: two peaks at \(\pm d/2\) give \(\cos(\pi x d/\lambda z)\) under an envelope, a lateral period of \(2\lambda z/d\) (29.15); the unapodized transmit beam only sets the envelope. A longer \(\lambda_x\) means a smaller lateral phase change per unit motion, so lateral estimates are noisier, alias later, and degrade linearly with depth.
Coherent summation of \(n\) echoes with independent noise raises signal power \(n^{2}\) and noise power \(n\). The compounded frame rate, which is the slow-time sampling rate of ultrafast Doppler, is PRF/\(n\), so \(v_{\mathrm{Nyq}}\) falls by \(n\). Motion during the \(n\) transmits adds a phase error, another reason to keep \(n\) small for fast flow.
29.11Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Doppler equation | \(f_D = 2vf_0\cos\theta/c\) | 6.49 kHz per m/s at 5 MHz, \(\theta\) = 0 |
| Angle error | \(\delta v/v \approx \tan\theta\,\delta\theta\) | keep \(\theta \le 60°\) |
| Nyquist velocity | \(v_{\mathrm{Nyq}} = c\,\mathrm{PRF}/(4f_0\cos\theta)\) | 0.308 m/s at 5 MHz, 4 kHz |
| Depth limit | \(z_{\max} = c/(2\,\mathrm{PRF})\) | 19.3 cm at 4 kHz |
| Range–velocity product | \(v_{\max}z_{\max} = c^{2}/(8f_0\cos\theta)\) | choose \(f_0\) for deep fast flow |
| Poiseuille spectrum | flat on \([0, f_{D,\max}]\), mean \(f_{D,\max}/2\) | volume flow from mean velocity |
| Geometric broadening | \(f_D\cos(\theta \mp \alpha)/\cos\theta\) | peak-velocity overestimate |
| STFT resolution | \(\Delta f \approx 2\,\mathrm{PRF}/N_w\), \(\Delta t = N_w/\mathrm{PRF}\) | spectrogram window choice |
| Kasai | \(\hat v = \frac{c\,\mathrm{PRF}}{4\pi f_0}\arg\hat R(1)\) | colour flow velocity |
| Loupas | divide by \(f_{\mathrm{dem}} + \frac{f_s}{2\pi}\arg\hat R(1,0)\) | removes RF downshift bias |
| Spectral width | \(\hat\sigma_\omega^{2} = \frac{2}{T^{2}}(1 - \lvert\hat R(1)\rvert/\hat R(0))\) | variance (turbulence) map; noise-biased |
| Regression wall filter | \(f_{-3\,\mathrm{dB}} \approx (0.44 + 0.41p)\,\mathrm{PRF}/N\) | colour-flow clutter rejection |
| SVD clutter filter | \(\mathbf{S}_f = \mathbf{S}(\mathbf{I} - \mathbf{V}_K\mathbf{V}_K^{H})\) | ultrafast Doppler, fUS, ULM |
| Power Doppler CNR | \((P_b/P_n)\sqrt{N - K}\) | ensemble length design |
| Multi-angle vector | \(\hat{\mathbf{v}} = (\mathbf{A}^{T}\mathbf{A})^{-1}\mathbf{A}^{T}\hat{\mathbf{f}}\); ratio \(\cot\gamma\) | vector flow, angle-free velocity |
| Transverse oscillation | \(\lambda_x = 2\lambda z/d\) | lateral velocity from one beam pair |
| Dual PRF | \(v_{\mathrm{ext}} = \lambda/(4\lvert T_2 - T_1\rvert)\) | anti-aliasing |
29.12Interview and qualifying-exam questions
I1Where does the factor 2 in the Doppler equation come from, and how large is the shift for carotid flow? · screen
The scatterer is a moving receiver on the way in and a moving source on the way back, so the path changes at \(2v\cos\theta\) and the echo phase rotates at \(f_D = 2vf_0\cos\theta/c\). For 1 m/s peak systolic flow at 60° and 5 MHz: 3.25 kHz. Follow-up: "is it a frequency shift of the echo?" In PW it is a phase change from pulse to pulse, sampled at the PRF; the shift of the RF spectrum within one pulse (\(v/c \approx 10^{-3}\)) is far below the bandwidth and is not what is measured.
I2Why does PW Doppler alias and CW not? Give the maximum unaliased velocity for a carotid at 3 cm, 5 MHz, 60°. · screen
PW samples the slow-time signal at the PRF; CW demodulates continuously. \(\mathrm{PRF}_{\max} = c/(2z)\) = 25.7 kHz at 3 cm, so \(v_{\max} = c^{2}/(8f_0\cos\theta\,z)\) = 3.95 m/s (29.4): no problem for a normal carotid, but a severe stenosis jet can exceed it. Follow-up: "what if the gate were at 15 cm?" 0.79 m/s; use baseline shift, lower \(f_0\), high-PRF or CW.
I3Sketch the Doppler spectrum of laminar flow in a vessel and explain why the maximum-velocity envelope overestimates. · onsite
Uniform insonation of a parabolic profile gives a flat spectrum from 0 to \(f_{D,\max}\) (equal area per velocity bin); a thin beam or a central gate tilts it toward the maximum. The aperture spreads each velocity over \(f_D[\cos(\theta + \alpha), \cos(\theta - \alpha)]/\cos\theta\): at F# 2 and 60° the edge reaches 1.39 \(f_{D,\max}\), so an envelope tracer reads 20–40 % high. The mean is unbiased to first order. Follow-up: "how would you correct it?" Calibrate the envelope for the aperture, use a smaller receive aperture, or use the mean.
I4Derive the Kasai estimator. Why is it the mean frequency, and when is it biased? · onsite
\(R(T) = E[z_{n+1}z_n^{*}] = \int S(f)e^{j2\pi fT}df\); for a spectrum symmetric about \(\bar f\), \(R(T) = e^{j2\pi\bar fT}\int S(f')\cos(2\pi f'T)\,df'\), so \(\arg R(T) = 2\pi\bar fT\) (29.8). Biases: asymmetric spectra (small), clutter residue (pulls toward tissue), aliasing, and the RF downshift: using \(f_{\mathrm{dem}}\) instead of the mean RF frequency scales \(\hat v\) by \(\bar f_{\mathrm{RF}}/f_{\mathrm{dem}}\), −10 % for 4.5 vs 5 MHz. Follow-up: Loupas 2D autocorrelation (29.9), which also lowers the variance by normalizing each estimate by the local mean RF frequency of the same block (depth averaging of \(\hat R(1)\) helps both estimators).
I5The variance map lights up in a vein with laminar flow. What is going on? · onsite
Low SNR: noise is white, adds to \(R(0)\) only and reduces \(\lvert R(1)\rvert/R(0)\) by \(\mathrm{SNR}/(1+\mathrm{SNR})\) (29.10). At PRF 4 kHz and 10 dB SNR a 150 Hz spectrum reads 306 Hz. Slow venous flow is also close to the wall filter, which shapes the residual spectrum. Fixes: subtract a noise estimate from \(R(0)\), threshold the variance by power, average over more samples. Follow-up: geometric broadening and shear also widen the spectrum without turbulence.
I6Design a wall filter: N = 10, PRF 3 kHz, 5 MHz, tissue moving at 5 mm/s, clutter 40 dB above blood. · onsite
Tissue Doppler is 32.5 Hz. Regression filters attenuate it by 14 dB (\(p\) = 0), 36 dB (\(p\) = 1) and 61 dB (\(p\) = 2), so \(p\) = 2 is needed to bring the residue below blood. Its cutoff is 380 Hz (rule: 378 Hz), i.e. flow slower than 0.117 m/s along a vessel at 60° is lost. To recover slow flow: longer \(N\) (cutoff ∝ PRF/N), or a lower PRF if the Nyquist limit allows, or an adaptive (SVD, eigen) filter. Follow-up: "why not an IIR filter?" Its transient eats most of a 10-sample ensemble unless it is initialized (projection initialization).
I7Why does an SVD filter beat a high-pass filter for microvascular imaging, and how do you pick K? · research
Slow microvascular flow and tissue motion overlap in Doppler frequency (1 mm/s at 15 MHz is 19.5 Hz, breathing and pulsatility give tens of Hz), so no frequency cutoff separates them. They differ in space: tissue moves coherently and is low rank in the Casorati matrix, blood is not (29.12). A short-ensemble high-pass also leaks motion through the window sidelobes (m29_ultrafast_power_doppler.m: background 11.5 dB above noise). K: from the knee of the singular spectrum, from the spatial similarity of singular vectors (Baranger et al.), or by a power threshold; too small leaves flash, too large removes slow blood that has not decorrelated. Follow-up: large fields with non-uniform motion need block-wise SVD; the noise subspace can also be cut at the high end.
I8Derive how power Doppler sensitivity scales with ensemble length. What stops you from making N arbitrarily large? · research
The background value is a mean of \(N - K\) exponential variables: mean \(P_n\), relative spread \(1/\sqrt{N-K}\); hence \(\mathrm{CNR} \approx (P_b/P_n)\sqrt{N-K}\) (29.13), ×2 per ×4 in \(N\). Limits: the stationarity of the scene (motion moves vessels between pixels and raises the clutter rank), the temporal resolution needed (fUS wants a power value every 0.2–1 s), data rate and memory (\(P\times N\) complex samples per block) and SVD cost (\(\propto PN^{2}\)). Follow-up: compounding adds \(10\log_{10}n_a\) dB to \(P_b/P_n\) at the cost of frame rate.
I9How does cross-beam vector Doppler work, and why is the lateral component noisier? Compare with transverse oscillation. · onsite
Each beam pair measures a projection \(f_k\) (29.14); least squares gives \(\mathbf{v}\) with covariance \(\sigma_f^{2}(\mathbf{A}^{T}\mathbf{A})^{-1}\). For \(\pm\gamma\) the lateral row coefficient is \(\propto\sin\gamma\), the axial \(\propto\cos\gamma\): ratio \(\cot\gamma\), 5.7 at 10°. Transverse oscillation creates a lateral wavelength \(\lambda_x = 2\lambda z/d\) (1.85 mm at 30 mm with \(d\) = 10 mm, six times \(\lambda\)), so the lateral phase per unit motion is six times smaller, the same geometric penalty in another form. Follow-up: both degrade with depth because the usable angle shrinks as \(d/z\).
I10How do sound-speed errors and aberration enter Doppler, and how is Doppler used in ULM? · research
\(v = cf_D/(2f_0\cos\theta)\) scales with the assumed \(c\) (a 3 % error is a 3 % velocity error), refraction changes the true \(\theta\), and defocusing (wrong \(c\), aberrating layers) widens the beam, so more angles contribute (more broadening) and the compounded SNR drops (Modules 32, 33). Phase aberration that is constant over the ensemble cancels in \(\arg R(1)\), which is why Doppler is robust to it while B-mode is not. In ULM the SVD filter of 29.6 isolates bubbles; the sign of the bubble's Doppler (or the direction of its track) separates up- and down-going flow, and track velocities give the velocity maps. Follow-up: motion between compounded angles creates a velocity-dependent phase error that biases both Doppler and ULM localization.
29.13MATLAB lab
Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run.
m29_doppler_spectrum_profile.mDoppler spectrum of plug, parabolic and blunted profiles; gate and elevation beam; geometric broadening; simulated periodogram; STFT resolution · change θ, vmax, F#, gate, Nw · runtime ≈ 4 s
%% m29_doppler_spectrum_profile.m
% Purpose : Doppler power spectrum of flow in a straight vessel. The spectrum
% is the histogram of f_D = 2 v f0 cos(theta)/c over the scatterers
% in the sample volume, broadened by the spread of beam angles
% across the aperture (geometric = transit-time broadening).
% Part 1: plug, parabolic and blunted profiles, full-lumen gate.
% Part 2: gate size and elevation beam (thin vs wide).
% Part 3: geometric broadening and the max-velocity overestimate.
% Part 4: simulated slow-time signal + averaged periodogram check.
% Part 5: STFT window length vs frequency/time resolution.
% Change : theta, vmax, a, F-number (alpha), gate length, Nw.
% Expect : f0 = 5 MHz, theta = 60 deg, vmax = 0.8 m/s: fDmax = 2597 Hz.
% Parabolic profile, whole lumen, wide elevation beam: flat
% spectrum, mean = fDmax/2 = 1299 Hz. Thin elevation beam: the
% spectrum rises toward fDmax. F# 2 (alpha = 14.0 deg): the
% broadened spectrum extends to fDmax*(1 + alpha*tan(theta)) =
% 1.42 fDmax; at -20 dB of its plateau it ends near 1.33 fDmax.
% Hann window, PRF 5 kHz: Nw = 128 -> bin 39.1 Hz, -6 dB width
% 2.00 bins = 78 Hz, 25.6 ms.
rng(1);
%% Parameters
c = 1540; % speed of sound (m/s)
f0 = 5e6; % centre frequency (Hz)
theta = 60*pi/180; % beam-to-flow angle (rad)
a = 3e-3; % vessel radius (m)
vmax = 0.8; % centreline velocity (m/s)
Fnum = 2; % F-number of tx and rx apertures
alpha = atan(1/(2*Fnum)); % aperture half-angle (rad)
fDmax = 2*vmax*f0*cos(theta)/c;
fprintf('fDmax = %.0f Hz (v = %.2f m/s, f0 = %.1f MHz, theta = %d deg)\n', fDmax, vmax, f0/1e6, round(theta*180/pi));
nb = 400; % frequency bins over [-0.1, 1.6] fDmax
fax = linspace(-0.1, 1.6, nb) * fDmax; df = fax(2) - fax(1);
hist1 = @(f, w) accumarray(min(max(round((f(:) - fax(1))/df) + 1, 1), nb), w(:), [nb 1])' / (sum(w(:)) * df);
% Cross-section coordinates: u = in-plane distance from the axis (u = s sin(theta),
% s along the beam), y = elevation. Uniform fine grid over the square, masked to the lumen.
ng = 801; [U, Y] = meshgrid(linspace(-a, a, ng)); R2 = (U.^2 + Y.^2) / a^2; in = R2 < 1;
prof = {@(r2) ones(size(r2)), @(r2) 1 - r2, @(r2) 1 - r2.^2}; % plug, parabolic, blunted (n = 4)
pname = {'plug', 'parabolic', 'blunted n=4'};
%% Part 1: profile shape, full-lumen gate, wide elevation beam
S1 = zeros(3, nb);
fprintf('Part 1 (whole lumen, wide elevation beam):\n');
for k = 1:3
v = vmax * prof{k}(R2(in)); f = 2*v*f0*cos(theta)/c;
S1(k, :) = hist1(f, ones(size(f)));
fprintf(' %-12s mean f = %6.0f Hz (%.3f fDmax) mean v = %.3f m/s\n', pname{k}, mean(f), mean(f)/fDmax, mean(v));
end
%% Part 2: gate length and elevation beam (parabolic profile)
fprintf('Part 2 (parabolic):\n');
gates = [1 0.3]; S2 = zeros(4, nb); ii = 0;
for elev = [1 0] % 1 = wide elevation beam, 0 = thin (y = 0 only)
for g = gates % gate covers |u| < g*a (centred)
ii = ii + 1;
if elev, m = in & abs(U) < g*a; uu = U(m); yy = Y(m);
else, uu = linspace(-g*a, g*a, 4001)'; yy = zeros(size(uu)); end
f = 2*vmax*(1 - (uu.^2 + yy.^2)/a^2)*f0*cos(theta)/c;
S2(ii, :) = hist1(f, ones(size(f)));
enames = {'thin', 'wide'};
fprintf(' elevation %-5s gate %3.0f %% of lumen: mean f = %.3f fDmax, min f = %.3f fDmax\n', ...
enames{elev + 1}, 100*g, mean(f)/fDmax, min(f)/fDmax);
end
end
%% Part 3: geometric broadening (tx and rx angles uniform in theta +/- alpha)
% f = v f0 (cos(theta_t) + cos(theta_r))/c ; linearized: triangle of half-width alpha*tan(theta)*f_D
v = vmax * prof{2}(R2(in)); v = v(randperm(numel(v), 20000)); % random subset of the lumen
fD = 2*v*f0*cos(theta)/c;
nd = 41; d = linspace(-alpha, alpha, nd); [Dt, Dr] = meshgrid(d);
fac = (cos(theta + Dt(:)) + cos(theta + Dr(:))) / (2*cos(theta)); % exact angle factor
fB = fD(:) * fac'; % every scatterer x angle pair
S3 = hist1(fB(:), ones(numel(fB), 1));
plateau = median(S1(2, fax > 0.2*fDmax & fax < 0.8*fDmax));
iEnd = find(S3 > 0.01*plateau, 1, 'last');
fprintf('Part 3: alpha = %.1f deg (F# %d): full extent = %.3f fDmax (exact cos(theta-alpha)/cos(theta) = %.3f, linear 1 + alpha tan(theta) = %.3f)\n', ...
alpha*180/pi, Fnum, max(fB(:))/fDmax, cos(theta - alpha)/cos(theta), 1 + alpha*tan(theta));
fprintf(' -20 dB envelope edge = %.3f fDmax -> apparent vmax = %.3f m/s (true %.2f)\n', fax(iEnd)/fDmax, vmax*fax(iEnd)/fDmax, vmax);
fprintf(' mean f unchanged? %.3f vs %.3f fDmax\n', mean(fB(:))/fDmax, mean(fD)/fDmax);
%% Part 4: simulated slow-time signal and averaged periodogram (parabolic, broadened)
% Each segment is an independent realization of K scatterers (uniform over the disk,
% random tx/rx angles within +/- alpha, complex Gaussian amplitudes) plus white noise.
PRF = 8e3; Nw = 256; nseg = 200; K = 300; t = (0:Nw-1)/PRF;
win = 0.5 - 0.5*cos(2*pi*(0:Nw-1)'/Nw); P = zeros(Nw, 1); fsum = 0;
for s = 1:nseg
rr = sqrt(rand(K,1));
fk = 2*vmax*(1 - rr.^2)*f0*cos(theta)/c .* ((cos(theta + alpha*(2*rand(K,1)-1)) + cos(theta + alpha*(2*rand(K,1)-1))) / (2*cos(theta)));
amp = (randn(1,K) + 1j*randn(1,K))/sqrt(2);
seg = (amp * exp(2j*pi*fk*t)).'; % Nw x 1 slow-time IQ
seg = seg + sqrt(K*1e-3) * (randn(Nw,1) + 1j*randn(Nw,1))/sqrt(2); % SNR 30 dB
P = P + abs(fft(seg .* win)).^2; fsum = fsum + mean(fk);
end
P = fftshift(P); fP = (-Nw/2:Nw/2-1)'*PRF/Nw;
fprintf('Part 4: periodogram mean f = %.0f Hz; model mean = %.0f Hz (K = %d scatterers x %d segments)\n', ...
sum(fP.*P)/sum(P), mean(fB(:)), K, nseg);
%% Part 5: STFT resolution (Hann window)
fprintf('Part 5: STFT at PRF = 5 kHz, Hann window\n Nw bin(Hz) -6dB width(Hz) dt(ms) dv at theta=0 (mm/s)\n');
PRF5 = 5e3;
for Nw = [32 64 128 256]
w = 0.5 - 0.5*cos(2*pi*(0:Nw-1)'/Nw); Z = 64*Nw;
W = abs(fft(w, Z)); W = W / W(1); k6 = find(W < 0.5, 1) - 1; % first bin below -6 dB
w6 = 2 * k6 / Z * PRF5; % two-sided -6 dB width (Hz)
fprintf(' %4d %6.1f %6.1f %5.1f %6.1f\n', Nw, PRF5/Nw, w6, Nw/PRF5*1e3, c*w6/(2*f0)*1e3);
end
%% Figures
figure('Position', [50 50 1100 380]);
subplot(1,3,1); plot(fax/fDmax, S1*fDmax, 'LineWidth', 1.5); xlabel('f / f_{D,max}'); ylabel('normalized power density');
title('Profile shape (whole lumen)'); legend(pname, 'Location', 'north'); ylim([0 5]); xlim([-0.1 1.6]); grid on;
subplot(1,3,2); plot(fax/fDmax, S2*fDmax, 'LineWidth', 1.5); xlabel('f / f_{D,max}'); ylim([0 8]); xlim([-0.1 1.6]); grid on;
title('Parabolic: gate and elevation beam'); legend('wide, 100%', 'wide, 30%', 'thin, 100%', 'thin, 30%', 'Location', 'northwest');
subplot(1,3,3); plot(fax/fDmax, S1(2,:)*fDmax, fax/fDmax, S3*fDmax, 'LineWidth', 1.5); hold on;
pl = median(P(fP > 0.2*fDmax & fP < 0.6*fDmax)); pm = median(S3(fax > 0.2*fDmax & fax < 0.6*fDmax));
plot(fP/fDmax, P/pl*pm*fDmax, 'k.-'); xlabel('f / f_{D,max}'); grid on; xlim([-0.1 1.6]);
title(sprintf('Broadening, F#%d, \\theta = %d deg', Fnum, round(theta*180/pi))); legend('no broadening', 'model', 'simulated periodogram');
What to look for: \(f_{D,\max}\) = 2597 Hz; mean frequency 1.000, 0.500 and 0.667 \(f_{D,\max}\) for plug, parabolic and blunted; with a thin beam and a 30 % gate the mean is 0.970 and the minimum 0.910; at F# 2 the broadened extent is 1.390 \(f_{D,\max}\) (linear 1.424), the −20 dB edge 1.336, an apparent \(v_{\max}\) of 1.07 m/s for a true 0.80; Hann STFT at PRF 5 kHz: 79.3 Hz −6 dB width for \(N_w\) = 128 (25.6 ms).
m29_kasai_loupas.mKasai versus Loupas on moving speckle with an RF downshift; variance estimate and its noise bias; aliasing curve · change fc, N, M, SNR, σv · runtime ≈ 7 s
%% m29_kasai_loupas.m
% Purpose : Kasai (1D, lag-one slow-time) versus Loupas (2D) autocorrelation
% velocity estimators on simulated IQ ensembles of moving speckle.
% The received pulse centre frequency fc is lower than the
% demodulation frequency fdem (attenuation downshift). Kasai uses
% fdem in the Doppler equation and is biased by fc/fdem; Loupas
% estimates the mean RF frequency from the lag-one correlation along
% depth and removes the bias. Also: the variance (spectral width)
% estimate sigma^2 = (2/T^2)(1 - |R(1)|/R(0)) and its noise bias.
% Change : fc, N (ensemble), M (range samples), snrdB, sigv (velocity spread).
% Expect : fdem = 5 MHz, fc = 4.5 MHz, PRF = 4 kHz: v_Nyq = 0.308 m/s.
% For v = 0.15 m/s Kasai reads about 0.9 x v (-10 %), Loupas about
% 1.00 x v. Std falls as N grows (roughly 1/sqrt(N) with velocity
% spread). Noise inflates the variance estimate: with SNR s the
% ratio |R(1)|/R(0) is multiplied by s/(1+s).
rng(1);
%% Parameters
c = 1540; % speed of sound (m/s)
fdem = 5e6; % demodulation (nominal transmit) frequency (Hz)
fc = 4.5e6; % received centre frequency after attenuation (Hz)
B = 0.6; % -6 dB fractional bandwidth of the received pulse
fs = 20e6; % IQ sampling rate along depth (Hz)
PRF = 4e3; % pulse repetition frequency (Hz)
T = 1/PRF;
vNyq = c*PRF/(4*fdem);
fprintf('v_Nyq (with fdem) = %.3f m/s, lambda = %.3f mm\n', vNyq, c/fdem*1e3);
% Gaussian envelope with -6 dB bandwidth B*fc: |G(f)| = exp(-(f/sf)^2/2), sf from -6 dB half width
sf = B*fc/2 / sqrt(2*log(2)); % spectral std (Hz) of amplitude spectrum
st = 1/(2*pi*sf); % envelope std in time (s)
% IQ sample of scatterers at delays tau (K x 1) with amplitudes a, observed at fast time t (1 x M):
iq = @(t, tau, a) a.' * (exp(-((t - tau).^2)/(2*st^2)) .* exp(2j*pi*(fc - fdem)*(t - tau)) .* exp(-2j*pi*fdem*tau));
% (note: t - tau uses implicit expansion, K x M)
kasai = @(Z) c*PRF/(4*pi*fdem) * angle(sum(sum(Z(:, 2:end) .* conj(Z(:, 1:end-1)))));
loupas = @(Z) c/2 * (PRF/(2*pi) * angle(sum(sum(Z(:, 2:end) .* conj(Z(:, 1:end-1)))))) / ...
(fdem + fs/(2*pi) * angle(sum(sum(Z(2:end, :) .* conj(Z(1:end-1, :))))));
%% Simulation loop
z0 = 30e-3; t0 = 2*z0/c; % gate centre depth and its delay
dens = 20e3; % scatterers per metre (20 per mm)
nT = 400; % trials per setting
cases = [4 8 16]; v = 0.15; sigv = 0.02; M = 8; snrdB = 20;
fprintf('v = %.3f m/s, sigma_v = %.3f m/s, M = %d range samples, SNR = %d dB, %d trials\n', v, sigv, M, snrdB, nT);
fprintf(' N Kasai mean bias std Loupas mean bias std (m/s, %%)\n');
res = zeros(numel(cases), 4);
for ic = 1:numel(cases)
N = cases(ic); vk = zeros(nT,1); vl = vk;
t = t0 + ((0:M-1) - (M-1)/2)/fs; % fast-time samples (1 x M)
Lw = (M/fs + 8*st) * c/2 + N*T*(abs(v) + 3*sigv); % spatial window (m)
for it = 1:nT
K = max(1, round(dens*Lw));
x0 = z0 + (rand(K,1) - 0.5)*Lw; % initial depths (m)
vk_s = v + sigv*randn(K,1); % scatterer velocities toward the probe
a = (randn(K,1) + 1j*randn(K,1))/sqrt(2);
Z = zeros(M, N);
for n = 1:N
tau = 2*(x0 - vk_s*(n-1)*T)/c; % moving toward the probe: delay shrinks
Z(:, n) = iq(t, tau, a).';
end
ps = mean(abs(Z(:)).^2);
Z = Z + sqrt(ps/10^(snrdB/10)) * (randn(M,N) + 1j*randn(M,N))/sqrt(2);
vk(it) = kasai(Z); vl(it) = loupas(Z);
end
res(ic, :) = [mean(vk) std(vk) mean(vl) std(vl)];
fprintf(' %2d %.4f %+5.1f%% %.4f %.4f %+5.1f%% %.4f\n', N, mean(vk), 100*(mean(vk)/v - 1), std(vk), ...
mean(vl), 100*(mean(vl)/v - 1), std(vl));
end
fprintf('Expected Kasai scale fc/fdem = %.3f (bias %+.1f %%)\n', fc/fdem, 100*(fc/fdem - 1));
%% Variance estimator: pure slow-time Gaussian-spectrum signal, spectral std sigf, with noise
sigf = 150; Nv = 16; nV = 4000; m = (0:Nv-1)';
fprintf('Variance estimate (Gaussian spectrum, sigma_f = %d Hz, N = %d):\n SNR(dB) sigma_f est (Hz) predicted with noise bias\n', sigf, Nv);
% Correlated Gaussian process by spectral shaping of white noise (long FFT, then crop)
Lf = 4096; ff = ((0:Lf-1) - Lf/2)'/Lf*PRF; H = exp(-ff.^2/(4*sigf^2)); H = ifftshift(H);
for snrdB2 = [40 20 10 0]
s2 = zeros(nV,1); snr2 = 10^(snrdB2/10);
for it = 1:nV/40
W = (randn(Lf,40) + 1j*randn(Lf,40))/sqrt(2);
X = ifft(bsxfun(@times, fft(W), H)); X = X(1:Nv, :); X = X / sqrt(mean(abs(X(:)).^2));
X = X + (randn(Nv,40) + 1j*randn(Nv,40))/sqrt(2*snr2);
for j = 1:40
z = X(:, j); R1 = sum(z(2:end).*conj(z(1:end-1)));
R0 = (sum(abs(z(2:end)).^2) + sum(abs(z(1:end-1)).^2))/2; % symmetric lag-0 (less finite-N bias)
s2((it-1)*40 + j) = 2/T^2 * (1 - abs(R1)/R0);
end
end
rho = snr2/(1 + snr2) * exp(-(2*pi*sigf*T)^2/2);
fprintf(' %3d %6.1f %6.1f\n', snrdB2, sqrt(max(mean(s2),0))/(2*pi), sqrt(2/T^2*(1 - rho))/(2*pi));
end
%% Aliasing: Kasai reading across velocities
vv = linspace(-2.5*vNyq, 2.5*vNyq, 121); vr = zeros(size(vv)); N = 8;
t = t0 + ((0:M-1) - (M-1)/2)/fs;
for i = 1:numel(vv)
K = round(dens*(M/fs*c/2 + 8*st*c/2 + 0.002)); x0 = z0 + (rand(K,1)-0.5)*(M/fs*c/2 + 8*st*c/2 + 0.002);
a = (randn(K,1) + 1j*randn(K,1))/sqrt(2); Z = zeros(M, N);
for n = 1:N, Z(:, n) = iq(t, 2*(x0 - vv(i)*(n-1)*T)/c, a).'; end
vr(i) = kasai(Z);
end
%% Figures
figure('Position', [50 50 1000 360]);
subplot(1,2,1); errorbar(cases - 0.15, res(:,1), res(:,2), 'o-'); hold on;
errorbar(cases + 0.15, res(:,3), res(:,4), 's-'); plot([2 18], [v v], 'k--');
xlabel('ensemble length N'); ylabel('velocity estimate (m/s)'); legend('Kasai (f_{dem})', 'Loupas 2D', 'true v', 'Location', 'southeast');
title(sprintf('f_c = %.1f MHz, f_{dem} = %.1f MHz', fc/1e6, fdem/1e6)); grid on; xlim([2 18]);
subplot(1,2,2); plot(vv, vr, '.', vv, vv*fc/fdem, 'k:'); xlabel('true velocity (m/s)'); ylabel('Kasai estimate (m/s)');
title('Aliasing: estimates wrap at \pm v_{Nyq}'); grid on;
What to look for: for \(v\) = 0.15 m/s Kasai reads 0.134–0.136 m/s (−9.6 to −10.4 %, expected −10.0 %) and Loupas 0.150 m/s (within 0.3 %); the standard deviation falls from 0.018 to 0.012 m/s (Kasai) and 0.014 to 0.010 m/s (Loupas) between \(N\) = 4 and 16; the width estimate of a 150 Hz spectrum reads 153, 187, 347 and 658 Hz at 40, 20, 10 and 0 dB SNR, against 149, 173, 306 and 645 Hz predicted.
m29_color_flow_wallfilter.mColour-flow frame: regression wall filters p = 0–3, Kasai velocity, power threshold, normalized variance; measured cutoffs against the rule · change p, N, PRF, CBR, vt · runtime ≈ 2 s
%% m29_color_flow_wallfilter.m
% Purpose : Colour-flow processing chain on a simulated 2D ensemble:
% clutter (40 dB above blood, slowly moving tissue) + blood in an
% inclined vessel (parabolic profile) + white noise ->
% polynomial-regression wall filter of degree p -> Kasai velocity,
% power and normalized variance (1 - |R(1)|/R(0)) per pixel ->
% power threshold. Also the filter's measured -3 dB cutoff,
% compared with (0.44 + 0.41 p) PRF/N from Module 09.
% Change : p (wall filter degree), N, PRF, CBR, tissue velocity vt.
% Expect : N = 12, PRF = 2 kHz, f0 = 5 MHz: v_Nyq = 0.154 m/s. p = 0 leaves
% tissue moving at 3 mm/s partly unfiltered (residual clutter >> blood,
% velocities pulled toward the tissue value, flash outside the vessel);
% p = 1-2 is the usable range; p = 3 eats slow flow near the walls.
rng(1);
%% Parameters
c = 1540; f0 = 5e6; PRF = 2e3; T = 1/PRF; N = 12;
theta = 60*pi/180; % beam (vertical) to vessel angle
a = 3e-3; vmax = 0.25; % vessel radius (m), centreline speed (m/s)
CBRdB = 40; SNRdB = 15; % clutter-to-blood and blood-to-noise power ratios (dB)
vt = 3e-3; % axial tissue velocity (m/s), toward the probe
sigb = 40; sigt = 3; % blood and tissue spectral std (Hz)
nz = 60; nx = 100; dx = 0.4e-3; % pixels (depth x lateral), pixel size (m)
vNyq = c*PRF/(4*f0);
fprintf('v_Nyq = %.3f m/s; axial centreline v = %.3f m/s; tissue Doppler = %.1f Hz\n', vNyq, vmax*cos(theta), 2*vt*f0/c);
[X, Z] = meshgrid(((1:nx) - nx/2)*dx, ((1:nz) - nz/2)*dx);
ang = pi/2 - theta; % vessel inclination from horizontal
dist = -X*sin(ang) + Z*cos(ang); % signed distance from the vessel axis
lumen = abs(dist) < a;
vax = zeros(nz, nx); vax(lumen) = vmax*(1 - (dist(lumen)/a).^2)*cos(theta); % axial speed toward probe
fb = 2*vax*f0/c; ft = 2*vt*f0/c;
P = nz*nx; n = (0:N-1);
% Signals (P pixels x N samples): sums of Kp phasors give Gaussian-like spectra
Kp = 8;
blood = zeros(P, N); tiss = zeros(P, N);
for k = 1:Kp
blood = blood + (randn(P,1) + 1j*randn(P,1))/sqrt(2*Kp) .* exp(2j*pi*(fb(:) + sigb*randn(P,1))*n*T);
tiss = tiss + (randn(P,1) + 1j*randn(P,1))/sqrt(2*Kp) .* exp(2j*pi*(ft + sigt*randn(P,1))*n*T);
end
blood(~lumen(:), :) = 0;
noise = (randn(P,N) + 1j*randn(P,N))/sqrt(2) * 10^(-SNRdB/20);
Xsig = 10^(CBRdB/20)*tiss + blood + noise;
%% Regression filters, responses and processing
fax = linspace(-PRF/2, PRF/2, 801); E = exp(2j*pi*n'*fax*T); % N x F test tones
fprintf(' p -3dB cutoff (Hz) approx(0.44+0.41p)PRF/N v_min(m/s) resid CBR(dB) detect%% flash%% mean|err| (m/s)\n');
maps = cell(1, 4);
for p = 0:3
V = bsxfun(@power, ((0:N-1)' - (N-1)/2)/N, 0:p); [Q, ~] = qr(V, 0); % orthonormal polynomial basis
A = eye(N) - Q*Q'; % regression (projection) filter
H2 = sum(abs(A*E).^2, 1)/N; % |H(f)|^2 as in Module 09
pos = fax >= 0; fp = fax(pos); hp = H2(pos);
fc3 = fp(find(hp >= 0.5, 1)); % first frequency reaching -3 dB
Y = Xsig * A.'; % filter every pixel (rows)
R0 = (sum(abs(Y(:,2:end)).^2, 2) + sum(abs(Y(:,1:end-1)).^2, 2))/2;
R1 = sum(Y(:,2:end).*conj(Y(:,1:end-1)), 2);
vk = c*PRF/(4*pi*f0)*angle(R1);
nfl = 10^(-SNRdB/10) * (N - p - 1)/N * (N-1); % expected noise-only R0
dmask = R0 > 2*nfl; % 3 dB above noise floor
clutRes = mean(sum(abs(10^(CBRdB/20)*tiss*A.').^2, 2)); bl = mean(sum(abs(blood(lumen(:),:)*A.').^2, 2));
inl = lumen(:) & dmask;
fprintf(' %d %6.0f %6.0f %.4f %6.1f %5.1f %5.1f %.4f\n', p, fc3, (0.44 + 0.41*p)*PRF/N, ...
c*fc3/(2*f0*cos(theta)), 10*log10(clutRes/bl), 100*mean(dmask(lumen(:))), 100*mean(dmask(~lumen(:))), mean(abs(vk(inl) - vax(inl))));
vm = vk; vm(~dmask) = -1.1*vNyq; maps{p+1} = reshape(vm, nz, nx); % below colour range: black
if p == 2, Hp2 = H2; end
end
fprintf('Detection threshold: R0 > 2 x expected noise power. Velocities in m/s are axial (toward probe > 0).\n');
fprintf('Slowest detectable flow (p = 2): about %.3f m/s along the vessel at theta = 60 deg\n', c*(0.44+0.82)*PRF/N/(2*f0*cos(theta)));
% Normalized variance map (p = 2): high where noise dominates
V = bsxfun(@power, ((0:N-1)' - (N-1)/2)/N, 0:2); [Q, ~] = qr(V, 0); A = eye(N) - Q*Q'; Y = Xsig*A.';
R0 = (sum(abs(Y(:,2:end)).^2, 2) + sum(abs(Y(:,1:end-1)).^2, 2))/2; R1 = sum(Y(:,2:end).*conj(Y(:,1:end-1)), 2);
nv = 1 - abs(R1)./R0;
fprintf('Normalized variance 1-|R1|/R0 (p = 2): lumen centre %.2f, lumen edge %.2f, outside (noise) %.2f\n', ...
median(nv(abs(dist(:)) < a/3)), median(nv(abs(dist(:)) > 0.75*a & lumen(:))), median(nv(~lumen(:))));
%% Figures
figure('Position', [40 40 1150 600]);
for p = 0:3
subplot(2,3,p+1); imagesc(X(1,:)*1e3, Z(:,1)*1e3, maps{p+1}, [-1.1*vNyq vNyq]); axis image; colorbar;
title(sprintf('Kasai v (m/s), regression p = %d', p)); xlabel('x (mm)'); ylabel('z (mm)');
end
colormap([0 0 0; jet(126)]); % first entry (black) = below threshold
subplot(2,3,5); plot(fax, 10*log10(Hp2 + 1e-12), 'LineWidth', 1.5); hold on; plot(ft*[1 1], [-80 5], 'k--');
xlabel('Doppler frequency (Hz)'); ylabel('|H(f)|^2 (dB)'); title('Regression p = 2, N = 12'); ylim([-80 5]); grid on;
subplot(2,3,6); imagesc(X(1,:)*1e3, Z(:,1)*1e3, reshape(nv, nz, nx), [-0.01 1]); axis image; colorbar;
title('1 - |R(1)|/R(0), p = 2'); xlabel('x (mm)'); ylabel('z (mm)');
What to look for: cutoffs 75, 145, 210, 278 Hz (rule 73, 142, 210, 278); residual clutter-to-blood 26.7, 6.9, −16.3, −42.0 dB for \(p\) = 0 … 3; flash outside the vessel 100, 99.2, 27.4, 1.8 %; mean velocity error 81, 58, 6.2 and 7.3 mm/s; normalized variance 0.10 at the lumen centre, 0.35 at the edge and 0.71 in noise.
m29_ultrafast_power_doppler.mSVD clutter filter of a simulated ultrafast ensemble, power Doppler CNR versus N against √(N − K), SVD versus FFT high-pass · change Am, CBRdB, BNRdB, fcut, vessel speeds · runtime ≈ 5 s
%% m29_ultrafast_power_doppler.m
% Purpose : Ultrafast power Doppler on a simulated slow-time ensemble.
% Clutter: speckle tissue that moves axially (3.3 Hz, 50 um, with a
% depth gradient so the motion is not rigid) at 40 dB above blood.
% Blood: five small vessels (150 um wide) with axial speeds from
% 1 to 30 mm/s, each pixel a complex Gaussian process with a
% Gaussian Doppler spectrum. White noise 6 dB above blood.
% Part 1: SVD clutter filter on the Casorati matrix (pixels x time),
% clutter rank Kc = singular powers > 10 dB above the median.
% Part 2: power Doppler CNR per vessel versus ensemble length N;
% theory: CNR = (P_b / P_n) sqrt(N - Kc) for independent noise.
% Part 3: SVD versus an ideal FFT high-pass (|f| > 75 Hz) at N = 200.
% Change : f0, PRF, motion amplitude Am, CBRdB, BNRdB, fcut, the vessel speeds.
% Expect : f0 = 15 MHz, PRF = 1 kHz, peak tissue Doppler 20 Hz (30 Hz at the
% bottom). Kc = 4 at N = 400; the first 4 components hold 99.9993 %
% of the tissue energy. CNR follows sqrt(N - Kc): 2.47 -> 5.00 from
% N = 100 to 400 (measured ratio 1.98). Slow vessels fall below theory
% at N = 25-50 (blood not yet decorrelated: it shares the tissue
% subspace). The FFT high-pass leaks tissue through the rectangular
% window (motion not periodic in the ensemble): background 11.5 dB
% above noise, no vessel detectable, while SVD leaves 0.99 x noise.
rng(1);
%% Parameters
c = 1540; f0 = 15e6; PRF = 1e3; T = 1/PRF; % compounded frame rate (Hz)
Nmax = 400; % longest ensemble (frames)
nz = 64; nx = 64; dpix = 50e-6; % pixels and pixel size (m)
P = nz*nx; t = (0:Nmax-1)*T;
CBRdB = 40; BNRdB = -6; % clutter/blood and blood/noise per sample (dB)
Am = 50e-6; fm = 3.3; grad = 0.5; % motion amplitude (m), rate (Hz, not periodic in the window), depth gradient
fcut = 75; % FFT high-pass cutoff for comparison (Hz)
vax = [1 3 10 20 30]*1e-3; % axial blood speeds (m/s)
zc = [8 20 32 44 56]; % vessel centre rows (pixels)
fD = 2*vax*f0/c;
fprintf('lambda = %.1f um; Doppler of the five vessels: %s Hz (Nyquist %d Hz)\n', c/f0*1e6, mat2str(round(fD*10)/10), PRF/2);
fprintf('Peak tissue velocity %.2f mm/s -> peak tissue Doppler %.1f Hz (x%.1f at the bottom)\n', 2*pi*fm*Am*1e3, 2*2*pi*fm*Am*f0/c, 1 + grad);
%% Signals
[XX, ZZ] = meshgrid(1:nx, 1:nz);
lab = zeros(nz, nx); % vessel label per pixel (0 = background)
for k = 1:5, lab(abs(ZZ - zc(k) - 0.15*(XX - nx/2)) <= 1.5) = k; end % slightly oblique, 3 px wide
% Tissue: fixed speckle amplitude x motion phase 4 pi f0 d(z,t)/c
depthw = 1 + grad*(ZZ(:) - 1)/(nz - 1); % P x 1 displacement weight
d = Am*sin(2*pi*fm*t + 0.7); % 1 x N displacement (m)
spk = (randn(P,1) + 1j*randn(P,1))/sqrt(2);
tiss = spk .* exp(1j*4*pi*f0/c * depthw * d); % P x N, unit power
% Blood: white noise shaped by a Gaussian spectrum centred at fD (per vessel)
Lf = 1024; ff = ((0:Lf-1) - Lf/2)/Lf*PRF; blood = zeros(P, Nmax);
for k = 1:5
idx = find(lab(:) == k); sb = 20 + 0.3*fD(k); % spectral std (Hz)
fa = mod(ff - fD(k) + PRF/2, PRF) - PRF/2; % aliased offset from fD
H = ifftshift(exp(-fa.^2/(4*sb^2))); % amplitude filter -> power std sb
W = (randn(numel(idx), Lf) + 1j*randn(numel(idx), Lf))/sqrt(2);
B = ifft(bsxfun(@times, fft(W, [], 2), H), [], 2);
B = B(:, 1:Nmax); blood(idx, :) = B / sqrt(mean(abs(B(:)).^2));
end
noise = (randn(P, Nmax) + 1j*randn(P, Nmax))/sqrt(2) * 10^(-BNRdB/20);
S = 10^(CBRdB/20)*tiss + blood + noise; % Casorati matrix, P x N
%% Part 1: SVD of the full ensemble
[U, Sg, V] = svd(S, 'econ'); sv = diag(Sg);
s2 = sv.^2; floorP = median(s2);
Kc = find(s2 > 10*floorP, 1, 'last'); % clutter: > 10 dB above the median singular power
fprintf('Part 1: N = %d, clutter rank Kc = %d (singular power > 10 dB above median); s1/s_median = %.1f dB\n', Nmax, Kc, 10*log10(s2(1)/floorP));
fprintf(' tissue energy in first Kc components: %.4f %% of total tissue energy left = %.2e\n', ...
100*(1 - norm(10^(CBRdB/20)*tiss - U(:,1:Kc)*(U(:,1:Kc)'*10^(CBRdB/20)*tiss), 'fro')^2/norm(10^(CBRdB/20)*tiss, 'fro')^2), ...
norm(10^(CBRdB/20)*tiss - U(:,1:Kc)*(U(:,1:Kc)'*10^(CBRdB/20)*tiss), 'fro')^2/norm(10^(CBRdB/20)*tiss, 'fro')^2);
%% Part 2: power Doppler CNR vs ensemble length
Ns = [25 50 100 200 400]; bg = lab(:) == 0; cnr = zeros(numel(Ns), 5); Kn = zeros(size(Ns));
fprintf('Part 2: CNR = (mean PD vessel - mean PD background) / std PD background\n');
fprintf(' N Kc V1(1mm/s) V2(3) V3(10) V4(20) V5(30) theory (Pb/Pn) sqrt(N - Kc)\n');
for i = 1:numel(Ns)
N = Ns(i); Sn = S(:, 1:N);
[Un, Sgn, Vn] = svd(Sn, 'econ'); q = diag(Sgn).^2; K = find(q > 10*median(q), 1, 'last'); Kn(i) = K;
F = Sn - Un(:,1:K)*(Sgn(1:K,1:K)*Vn(:,1:K)');
PD = mean(abs(F).^2, 2);
for k = 1:5, cnr(i,k) = (mean(PD(lab(:)==k)) - mean(PD(bg))) / std(PD(bg)); end
fprintf(' %3d %2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n', N, K, cnr(i,:), 10^(BNRdB/10)*sqrt(N - K));
if N == 200, PDsvd = PD; F200 = F; end
end
fprintf('CNR ratio N = 400 vs 100 (mean of vessels): %.2f (sqrt(4) = 2)\n', mean(cnr(5,:)./cnr(3,:)));
%% Part 3: SVD vs FFT high-pass at N = 200
N = 200; Sn = S(:, 1:N);
fk = ((0:N-1) - N*(0:N-1 >= N/2))*PRF/N; % FFT bin frequencies (Hz)
G = fft(Sn, [], 2); G(:, abs(fk) <= fcut) = 0; Fh = ifft(G, [], 2);
PDhp = mean(abs(Fh).^2, 2);
cnrH = zeros(1,5); cnrS = zeros(1,5);
for k = 1:5
cnrH(k) = (mean(PDhp(lab(:)==k)) - mean(PDhp(bg))) / std(PDhp(bg));
cnrS(k) = (mean(PDsvd(lab(:)==k)) - mean(PDsvd(bg))) / std(PDsvd(bg));
end
resH = mean(PDhp(bg)) / (10^(-BNRdB/10)); resS = mean(PDsvd(bg)) / (10^(-BNRdB/10));
fprintf('Part 3 (N = 200): background PD / noise power: SVD %.2f, high-pass %.1f (%.1f dB above noise)\n', resS, resH, 10*log10(resH));
fprintf(' vessel: 1 mm/s 3 mm/s 10 mm/s 20 mm/s 30 mm/s\n SVD CNR %6.2f %6.2f %6.2f %6.2f %6.2f\n HP CNR %6.2f %6.2f %6.2f %6.2f %6.2f\n', cnrS, cnrH);
%% Figures
figure('Position', [40 40 1150 640]);
subplot(2,3,1); plot(10*log10(s2/s2(1)), '.-'); hold on; plot([Kc Kc], [-80 0], 'k--'); grid on;
xlabel('singular index'); ylabel('s_k^2 / s_1^2 (dB)'); title(sprintf('Singular spectrum, N = %d, Kc = %d', Nmax, Kc)); xlim([1 Nmax]);
subplot(2,3,2); plot(t(1:200)*1e3, real(V(1:200,1)), t(1:200)*1e3, real(V(1:200,Kc+5)), 'LineWidth', 1.2); grid on;
xlabel('slow time (ms)'); ylabel('Re v_k'); legend('k = 1 (tissue)', sprintf('k = %d (blood/noise)', Kc+5)); title('Temporal singular vectors');
subplot(2,3,3); loglog(Ns, cnr, 'o-'); hold on; loglog(Ns, 10^(BNRdB/10)*sqrt(Ns), 'k--', 'LineWidth', 1.5); grid on;
xlabel('ensemble length N'); ylabel('power Doppler CNR'); title('Sensitivity vs ensemble length'); legend('1', '3', '10', '20', '30 mm/s', 'sqrt(N) P_b/P_n', 'Location', 'northwest');
subplot(2,3,4); imagesc((1:nx)*dpix*1e3, (1:nz)*dpix*1e3, reshape(10*log10(mean(abs(S(:,1:200)).^2,2)), nz, nx)); axis image; colorbar; colormap(gray);
title('Unfiltered mean power (dB)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,5); imagesc((1:nx)*dpix*1e3, (1:nz)*dpix*1e3, reshape(10*log10(PDsvd*10^(BNRdB/10)), nz, nx), [-2 3]); axis image; colorbar;
title('Power Doppler, SVD, N = 200 (dB re noise)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,6); imagesc((1:nx)*dpix*1e3, (1:nz)*dpix*1e3, reshape(10*log10(PDhp*10^(BNRdB/10)), nz, nx), [-2 20]); axis image; colorbar;
title(sprintf('PD, high-pass > %d Hz (dB re noise)', fcut)); xlabel('x (mm)'); ylabel('z (mm)');
What to look for: clutter rank \(K\) = 4 at \(N\) = 400 (first four components hold 99.9993 % of the tissue energy, \(\sigma_1^{2}\) 58.5 dB above the median); CNR of the 10 mm/s vessel 1.23, 1.54, 2.61, 3.44, 4.68 for \(N\) = 25, 50, 100, 200, 400 against 1.20, 1.74, 2.47, 3.52, 5.00; mean ratio 1.98 from 100 to 400 frames; the 1 mm/s vessel lags at small \(N\) (0.41 at 25); at \(N\) = 200 the SVD background is 0.99 × noise and the high-pass background 11.5 dB above noise, with no vessel detectable.
m29_vector_flow_multiangle.mMulti-angle Doppler vectors: Monte Carlo against the least-squares covariance, cot γ and cot(γ/2) laws, transverse-oscillation wavelength from the Fraunhofer field · change γ, K, mode, SNR, N · runtime ≈ 2 s
%% m29_vector_flow_multiangle.m
% Purpose : 2D velocity vectors from multi-angle plane-wave Doppler.
% Each transmit/receive angle pair (alpha_k, beta_k) measures the
% projection f_k = -(f0/c) [ (sin a + sin b) vx + (cos a + cos b) vz ]
% (x lateral, z depth downward). K >= 2 pairs give a least-squares
% vector. Part 1: Monte Carlo (Kasai on simulated ensembles) versus
% the covariance sigma_f^2 (A'A)^-1. Part 2: the cross-beam cot(gamma)
% law for the lateral/axial precision ratio. Part 3: the lateral
% wavelength of transverse oscillation, lambda_x = 2 lambda z / d,
% measured from the receive field of a two-peak apodization.
% Change : gamma (angle span), K, mode (rx = tx or rx = 0), snrdB, N.
% Expect : precision ratio sigma_vx/sigma_vz = cot(gamma) for 2 symmetric
% pairs with rx = tx: 5.67 at 10 deg, 2.75 at 20 deg. With rx fixed
% at 0 deg the ratio is cot(gamma/2), about twice as large.
% lambda_x: z = 30 mm, d = 10 mm, f0 = 5 MHz -> 1.85 mm.
rng(1);
%% Parameters
c = 1540; f0 = 5e6; PRFtx = 10e3; % transmit rate (Hz); per-angle PRF = PRFtx/K
N = 32; snrdB = 10; nT = 2000; % ensemble per angle, per-sample SNR, trials
v = [0.20, -0.10]; % true (vx, vz) in m/s; vz < 0 moves toward the probe
Arow = @(al, be) -(f0/c) * [sin(al) + sin(be), cos(al) + cos(be)];
%% Part 1 and 2: Monte Carlo vs covariance
fprintf('True v = (%.2f, %.2f) m/s, N = %d per angle, SNR = %d dB, per-angle PRF = PRFtx/K\n', v(1), v(2), N, snrdB);
fprintf(' mode K gamma sd_vx MC sd_vx th sd_vz MC sd_vz th ratio cot(g) cot(g/2) v_Nyq,axial(m/s)\n');
modes = {'rx=tx', 'rx=0'};
for im = 1:2
for K = [2 5]
for g = [5 10 20 30]
gam = g*pi/180; al = linspace(-gam, gam, K)'; be = al * (im == 1);
A = zeros(K, 2); for k = 1:K, A(k, :) = Arow(al(k), be(k)); end
PRF = PRFtx/K; f = A*v'; m = (0:N-1)'; est = zeros(nT, 2); fh = zeros(K, nT);
for k = 1:K
z = exp(2j*pi*f(k)*m/PRF) * exp(2j*pi*rand(1, nT)) + (randn(N, nT) + 1j*randn(N, nT))/sqrt(2*10^(snrdB/10));
fh(k, :) = PRF/(2*pi) * angle(sum(z(2:end, :).*conj(z(1:end-1, :)), 1));
end
est = (A \ fh)'; % least squares, one vector per trial
sf = std(fh(1, :) - f(1)); % per-angle Doppler std (Hz)
Cv = sf^2 * inv(A'*A);
fprintf(' %-6s %d %2d %.4f %.4f %.4f %.4f %5.2f %5.2f %5.2f %.3f\n', modes{im}, K, g, std(est(:,1)), sqrt(Cv(1,1)), ...
std(est(:,2)), sqrt(Cv(2,2)), std(est(:,1))/std(est(:,2)), cot(gam), cot(gam/2), c*PRF/(4*f0));
if im == 1 && K == 2 && g == 10, cloud = est; end
end
end
end
fprintf('Mean estimate (rx=tx, K = 2, 10 deg): (%.3f, %.3f) m/s\n', mean(cloud(:,1)), mean(cloud(:,2)));
%% Part 3: transverse oscillation lateral wavelength from a two-peak receive apodization
lam = c/f0; d = 10e-3; w = 2.5e-3; % peak separation and Gaussian peak width (m)
xi = linspace(-20e-3, 20e-3, 4096); dxi = xi(2) - xi(1);
apod = exp(-((xi - d/2)/w).^2) + exp(-((xi + d/2)/w).^2);
fprintf('Transverse oscillation, f0 = %.0f MHz, d = %.0f mm:\n z (mm) lambda_x measured (mm) 2*lambda*z/d (mm)\n', f0/1e6, d*1e3);
for zf = [20 30 40]*1e-3
x = linspace(-6e-3, 6e-3, 1201);
h = exp(-2j*pi*x'*xi/(lam*zf)) * apod.' * dxi; % Fraunhofer receive field at depth zf (focused)
hr = real(h); % field is real (symmetric apodization)
big = abs(x(1:end-1)') < 3*lam*zf/(pi*w); % inside 3x the Gaussian envelope radius (zeros of h are the crossings, so do not mask on |h|)
zc = find(big & (hr(1:end-1) > 0 & hr(2:end) <= 0 | hr(1:end-1) < 0 & hr(2:end) >= 0));
pz = mean(diff(x(zc))) * 2; % period from zero crossings
fprintf(' %4.0f %6.3f %6.3f\n', zf*1e3, pz*1e3, 2*lam*zf/d*1e3);
if abs(zf - 30e-3) < 1e-6, xTO = x; hTO = h; end
end
%% Figures
figure('Position', [50 50 1100 360]);
subplot(1,3,1); plot(cloud(:,1), cloud(:,2), '.', 'MarkerSize', 3); hold on; plot(v(1), v(2), 'k+', 'MarkerSize', 14, 'LineWidth', 2);
axis equal; grid on; xlabel('v_x (m/s)'); ylabel('v_z (m/s)'); title('LS estimates, K = 2, \pm10 deg, rx = tx');
gg = linspace(2, 40, 100)*pi/180;
subplot(1,3,2); semilogy(gg*180/pi, cot(gg), 'LineWidth', 1.5); hold on; semilogy(gg*180/pi, cot(gg/2), '--', 'LineWidth', 1.5);
grid on; xlabel('\gamma (deg)'); ylabel('\sigma_{vx} / \sigma_{vz}'); legend('rx = tx: cot \gamma', 'rx = 0: cot(\gamma/2)'); title('Lateral precision penalty');
subplot(1,3,3); plot(xTO*1e3, real(hTO)/max(abs(hTO)), xTO*1e3, abs(hTO)/max(abs(hTO)), '--', 'LineWidth', 1.3);
grid on; xlabel('lateral x (mm)'); ylabel('receive field (norm.)'); title('TO field, z = 30 mm, d = 10 mm');
What to look for: Monte Carlo and covariance standard deviations agree to the second digit; for \(K\) = 2 with steered receive the ratio is 11.34, 5.66, 2.82, 1.68 at 5, 10, 20, 30° (cot γ: 11.43, 5.67, 2.75, 1.73), and with receive at 0° 22.83, 11.41, 5.56, 3.65 (cot(γ/2)); five angles spread uniformly over the same span lower both standard deviations but raise the ratio to about 1.4 cot γ (8.03 at 10°), because the inner angles add little lateral information; \(\lambda_x\) = 1.231, 1.849, 2.466 mm at 20, 30, 40 mm (theory 1.232, 1.848, 2.464).
29.14Self-check
Q17.5 MHz, θ = 45°, measured \(f_D\) = 2 kHz. What is v?
\(v = cf_D/(2f_0\cos\theta)\) \(= 1540 \times 2000\,/\,(2 \times 7.5\times 10^{6} \times 0.707)\) = 0.29 m/s.
Q2A gate at 8 cm, 3.5 MHz, θ = 0. Maximum PRF and unaliased velocity?
\(\mathrm{PRF}_{\max} = c/(2z)\) = 9.63 kHz; \(v_{\max} = c^{2}/(8f_0 z)\) = 1540²/(8 × 3.5 × 10⁶ × 0.08) = 1.06 m/s.
Q3Plug versus parabolic flow, same centreline speed: which has the higher mean Doppler frequency under uniform insonation, and by how much?
Plug: mean = \(f_{D,\max}\); parabolic: \(f_{D,\max}/2\). Plug is twice as high (29.5).
Q4The lag-one phase is 120° at PRF 2 kHz and 5 MHz. Velocity? What if the true velocity were larger than \(v_{\mathrm{Nyq}}\)?
\(\hat v = c\,\mathrm{PRF}\,\varphi/(4\pi f_0)\) \(= 1540 \times 2000 \times (2\pi/3)\,/\,(4\pi \times 5\times 10^{6})\) = 0.103 m/s (\(v_{\mathrm{Nyq}}\) = 0.154 m/s). The same phase also fits \(0.103 + 2m \times 0.154\) m/s, e.g. 0.411 or −0.205 m/s.
Q5N = 16, PRF 4 kHz, regression order 1. Cutoff, and the slowest axial velocity kept at 5 MHz?
\((0.44 + 0.41)\times 4000/16\) = 213 Hz; \(v = cf/(2f_0)\) = 3.3 cm/s axial (29.11).
Q6Blood-to-noise is −10 dB per sample after compounding and K = 5. What ensemble gives a power Doppler CNR of 3?
\(0.1\sqrt{N - 5} = 3\) → \(N\) = 905 frames, about 0.9 s at 1 kHz (29.13). Compounding 9 angles instead of 1 would add up to 9.5 dB (\(P_b/P_n\) = 0.9) and bring it to \(N\) = 16, at a ninth of the frame rate for the same PRF.
Q7Cross-beam with ±15° (receive steered), axial precision 2 mm/s. Lateral precision?
\(\cot 15°\) = 3.73, so 7.5 mm/s; with receive fixed at 0°, \(\cot 7.5°\) = 7.6, i.e. 15.2 mm/s (29.14).
29.15Figures from the literature
Measured data behind the models of this module. Figure numbers checked against the author manuscripts or the publisher's page; links only, no images reproduced.
The transverse oscillation field, and what it shows in a carotid bifurcation
J. A. Jensen, S. I. Nikolov, A. C. H. Yu and D. Garcia, "Ultrasound vector flow imaging, Part I: Sequential systems," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63(11), 1704–1721 (2016). doi:10.1109/TUFFC.2016.2600763
Look at Fig. 4 (the two-peak receive apodization and the left and right TO point-spread functions: the lateral oscillation of (29.15)) and Fig. 5 (TO vector flow image of a carotid bifurcation with a vortex near peak systole). Why: the vortex is invisible to axial Doppler, and the TO field is the mechanism that makes it measurable from one aperture.
Cross-beam geometry and Doppler bandwidth versus beam geometry
Same review as above (Part I). doi:10.1109/TUFFC.2016.2600763
Look at Fig. 1 (cross-beam acquisition, the \(K\) = 2 case of (29.14)) and Fig. 7 (Doppler bandwidth versus beam geometry, and a spectrogram from a transversely insonated common carotid). Why: Fig. 7 is the measured version of the geometric broadening of 29.3: at 90° the mean Doppler is zero but the spectrum still has a width proportional to the transverse velocity, which is the idea behind bandwidth-based transverse velocity estimation.
Plane-wave vector flow in vivo
J. A. Jensen, S. I. Nikolov, A. C. H. Yu and D. Garcia, "Ultrasound vector flow imaging, Part II: Parallel systems," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63(11), 1722–1732 (2016). doi:10.1109/TUFFC.2016.2598180
Look at Fig. 7 (plane-wave vector flow of the jugular vein and carotid artery at three points of the cardiac cycle) and Fig. 8 (high-frame-rate vector flow of a neonatal ventricular septal defect). Why: with plane waves every pixel is sampled continuously, so the multi-angle estimator of Demo 29.6 runs over the whole field at kHz rates.
From conventional Doppler to functional ultrasound in the rat brain
Same review (Part II). doi:10.1109/TUFFC.2016.2598180
Look at Fig. 9 (rat cerebral blood flow with conventional velocity imaging, plane-wave power Doppler and functional ultrasound). Why: the same brain imaged with a short conventional ensemble and with a long ultrafast ensemble: the gain in visible small vessels is the \(\sqrt{N}\) law of (29.13) plus the SVD filter.
Functional ultrasound: principle and first applications
E. Macé, G. Montaldo, I. Cohen, M. Baulac, M. Fink and M. Tanter, "Functional ultrasound imaging of the brain," Nature Methods 8, 662–664 (2011). doi:10.1038/nmeth.1641
Look at Fig. 1 (principles of fUS in the rat brain: plane-wave acquisition and power Doppler of blood volume) and Fig. 2 (applications: whisker-evoked activation in cortex and thalamus, and the spread of epileptic seizures). Why: it turned ultrafast power Doppler sensitivity into a neuroimaging modality, and it is the acquisition on which ULM of the brain was later built.
29.16Go deeper
- J. A. Jensen, Estimation of Blood Velocities Using Ultrasound: A Signal Processing Approach, Cambridge University Press (1996). The standard text: Doppler signal model, spectral estimation, autocorrelation and cross-correlation estimators, clutter.
- C. Kasai, K. Namekawa, A. Koyano and R. Omoto, "Real-time two-dimensional blood flow imaging using an autocorrelation technique," IEEE Trans. Sonics Ultrason. SU-32, 458–464 (1985). The colour-flow estimator of 29.4.
- T. Loupas, J. T. Powers and 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 TUFFC 42(4), 672–688 (1995). doi:10.1109/58.393110
- V. L. Newhouse, E. S. Furgason, G. F. Johnson and D. A. Wolf, "The dependence of ultrasound Doppler bandwidth on beam geometry," IEEE Trans. Sonics Ultrason. SU-27, 50–59 (1980). Geometric and transit-time broadening as one effect.
- J. Bercoff, G. Montaldo, T. Loupas, D. Savery, F. Mézière, M. Fink and M. Tanter, "Ultrafast compound Doppler imaging: providing full blood flow characterization," IEEE TUFFC 58(1), 134–147 (2011). doi:10.1109/TUFFC.2011.1780
- C. Demené et al., "Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity," IEEE Trans. Med. Imaging 34(11), 2271–2285 (2015). doi:10.1109/TMI.2015.2428634 The SVD filter of 29.6.
- J. A. Jensen and P. Munk, "A new method for estimation of velocity vectors," IEEE TUFFC 45(3), 837–851 (1998). doi:10.1109/58.677749 Transverse oscillation.
- Lab cross-links: lag-one phase estimation and speckle tracking (Module 11), SVD (Module 06), regression filters (Module 09); next: elastography (Module 30) reuses the same estimators on tissue, and ULM (Module 31) starts from the SVD-filtered ensemble.