A Bug That Doesn't Crash Is Worse Than One That Does — Building a QUS Python Refresher

I recently rebuilt my quantitative ultrasound (QUS) signal processing knowledge from scratch, in Python, as a five-notebook refresher repo. Along the way I reproduced a bug worth writing up on its own, because it never raised an exception — it just quietly returned zero.

The setup

Quantitative ultrasound tries to turn raw RF echo data into numbers that describe tissue — not just an image a radiologist reads, but a spectral slope or a backscatter coefficient that correlates with scatterer size, concentration, or fibrosis. The foundational method here is the Lizzi-Feleppa framework from 1983: fit a line to the log-power spectrum within a transducer's usable bandwidth, per depth window, and read off the slope, intercept, and midband value.

To get this back into muscle memory, I built qus-python-refresher — five Jupyter notebooks that simulate an RF echo from point scatterers and walk through the full pipeline: pulse generation, envelope detection, moving-window spectral analysis, and Lizzi-Feleppa parameter estimation.

The bug

Point-scatterer echo simulation needs two separate time axes:

• a short axis just long enough to hold the transmit pulse (t_pulse)
• a long axis covering the full imaging depth (t), onto which delayed, scaled copies of the pulse get summed

Building the pulse on the wrong axis — reusing t instead of t_pulse — is an easy mistake. The Gaussian envelope still decays to numerical zero away from its center, so if you plot the "pulse" it looks completely normal. The array itself, though, is now the same length as the full echo instead of a few dozen samples.

The consequence shows up in the placement loop's bounds check:

if idx + len(pulse) < len(rf):
    rf[idx:idx+len(pulse)] += a * pulse

With a pulse the same length as rf, idx + len(pulse) >= len(rf) for every idx >= 0 — the condition can never be true. In a 200-scatterer test I ran to confirm this: 0 out of 200 scatterers get placed with the buggy pulse, versus 200 out of 200 with the correctly-scoped one. rf comes back as a normal, correctly-shaped, all-zero array. No exception. No warning. Just silence.

That's the failure mode worth internalizing: shape mismatches that happen to broadcast are more dangerous than the ones that crash loudly. A flat line in a B-mode image would be an easy bug to notice; buried a few layers deep in a larger pipeline, it isn't.

What's in the repo

• 01 — RF Pulse and Envelope. A single transmit pulse, Hilbert transform, envelope, phase, log compression. No scatterers yet, on purpose — isolates the RF-to-B-mode pipeline before the second time axis shows up.
• 02 — Scatterer Echo Simulation. The bug above, reproduced with real numbers. Also: diffuse (Rayleigh-regime) vs. specular scattering, and frequency-dependent attenuation.
• 03 — Moving-Window Spectral Analysis. scipy.signal.spectrogram's default window is Tukey(0.25), not rectangular — easy to assume otherwise if you never pass window= explicitly. This notebook compares window lengths (1/3/10 µs — a real axial-resolution-vs-frequency-resolution trade-off) and window types (rectangular vs. Hann, and the leakage difference between them).
• 04 — Lizzi-Feleppa Slope Estimation. Slope, intercept, midband fit, per depth window. The interesting result here: scaling echo amplitude by 5x shifts intercept and midband by very close to the theoretical 20·log10(5) = 13.98 dB — but only once you exclude depth windows that sit before the first scatterer, where both runs are noise-floor-dominated and the shift is 0 dB regardless of gain. Averaging over the full depth range without that exclusion gave a misleading ~9.9 dB, and figuring out why was a good reminder to sanity-check estimator behavior against a known-zero baseline before trusting a summary statistic. Slope itself, unlike intercept/midband, comes out essentially unchanged by the gain — which is exactly why reference phantom normalization exists.
• 05 — Attenuation and Spectral Downshift. Why an uncompensated spectral slope is depth-biased: attenuation removes high frequencies faster than low ones, so the local spectrum's shape — not just its amplitude — shifts with depth. Includes a spectral-centroid-vs-depth demo (the observable the spectral shift method is built on), without implementing the full reference-phantom-based compensation.

Shared numerics (pulse generation, scatterer field simulation, attenuation, spectral fitting) live in a small src/qus_refresher package rather than being copy-pasted across notebooks, so a fix in one place doesn't drift out of sync elsewhere. Each notebook is stored twice — a jupytext "percent format" .py as the reviewable source, and an executed .ipynb with real, embedded plot outputs.

There's also a standalone NOMENCLATURE.md — a working reference on RF/QUS terminology, formulas, and trade-offs, and a reading list starting from Oelze & Mamou's 2016 IEEE TUFFC review, through the original Lizzi et al. 1983 JASA paper, to Yao/Zagzebski's work on reference-phantom backscatter coefficient extraction. I also wrote it up as its own post — Ultrasound RF & QUS — Nomenclature and Concepts — if you'd rather read the reference on its own, independent of this repo's code.

Not in scope yet: reference phantom normalization / backscatter coefficient extraction — most of the background papers behind this repo are specifically about that method, and it's the natural next notebook rather than something I decided to skip.

On how this got built

I built this with Claude doing a lot of the actual scaffolding and notebook-writing in a pair-programming style — I described the physics and the pedagogical structure I wanted, reviewed and pushed back on the design (the shared src/ package split, the notebook boundaries, catching that misleading 9.9 dB average above), and every notebook was actually executed end-to-end rather than hand-assembled, so the numbers and plots in the repo are real outputs, not illustrative guesses. It's a genuinely useful way to work through this kind of refresher fast, as long as you're willing to interrogate the numbers it hands you back.

Try it

Repo: qus-python-refresher — public. environment.yml + pip install -e . gets you a working environment; each notebook runs standalone once the shared package is installed. Feedback, issues, and "actually your slope-invariance argument is wrong because..." are all welcome.