diff --git a/CHANGELOG.md b/CHANGELOG.md
index 69da0f5..4752b68 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,54 @@ All notable changes to this project are documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+
+- `Latency()` on `SimpleResampler` and `SimpleResamplerFloat32`, returning the
+  startup deficit in output samples so a caller can prime a real-time FIFO
+  before the first output. (#51)
+- Streaming documentation (latency, the `Flush` contract, per-channel instance
+  requirements) and a real-time FIFO example at `examples/streaming`. (#51)
+
+### Fixed
+
+- Polyphase phase-boundary coefficient interpolation used a wrapped neighbor,
+  degrading THD+N at ratios with active sub-phase interpolation. The committed
+  measurement shows an 86.26 dB improvement (wrapped -54.46 dB versus flat
+  -140.72 dB) at 44100 to 64000, with similar magnitude at other
+  active-interpolation ratios measured during the investigation but not
+  committed as tests; exact-rational ratios such as 44100 to 48000 were
+  unaffected. (#51)
+- Severe non-integer downsampling (beyond roughly 1:16) corrupted output with
+  repeated stale samples and grew internal history without bound. (#51)
+- `Flush` over-padded each filter stage by one zero, emitting about 2 phantom
+  samples; `Process` plus `Flush` now totals within
+  `[floor(n*ratio), ceil(n*ratio)+1]`. (#51)
+- At unity ratio (`inputRate == outputRate`), `Process` returned the caller's
+  own input slice; it now returns an owned buffer. (#51)
+- Cubic (`QualityQuick`) resampling computed its first output segments from a
+  fictional zero history and never emitted the final segments; output is now
+  aligned to real data, with the first output after 2 input samples. Its
+  `Process` also returned an aliased empty slice during priming; it now returns
+  an owned buffer. (#51)
+- NaN sample rates are now rejected by all constructors, covering both the
+  `NewEngine`/`NewEngineFloat32` engine path and the `New(config)` pipeline
+  path (`New`, `NewMultiChannel`, `NewStereo`, `NewSimple`, and the preset
+  helpers). (#51)
+- Half-band stage construction errors now propagate instead of silently
+  substituting a nearest-neighbor stub. (#51)
+- `GetLatency` now accounts for decimation and cubic stages. (#51)
+
+### Changed
+
+- `Flush` is now terminal: a second `Flush` returns an empty slice, and a
+  `Process` call after `Flush` starts a fresh stream instead of convolving
+  against leftover padding. (#51)
+- `QualityQuick` through `NewEngine` and `NewEngineFloat32` now uses cubic
+  interpolation (matching `New()` and the documented contract) instead of a
+  full FIR pipeline; latency drops accordingly. (#51)
+
 ## [1.4.0] - 2026-05-29
 
 ### Added
@@ -69,6 +117,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   design, quality presets, multi-channel and streaming support, validated against
   libsoxr.
 
+[Unreleased]: https://github.com/tphakala/go-audio-resampler/compare/v1.4.0...HEAD
 [1.4.0]: https://github.com/tphakala/go-audio-resampler/compare/v1.3.0...v1.4.0
 [1.3.0]: https://github.com/tphakala/go-audio-resampler/compare/v1.2.0...v1.3.0
 [1.2.0]: https://github.com/tphakala/go-audio-resampler/compare/v1.1.0...v1.2.0
diff --git a/README.md b/README.md
index f1ea461..210c7ad 100644
--- a/README.md
+++ b/README.md
@@ -112,11 +112,12 @@ import (
 )
 
 func main() {
-    // Create a resampler for CD to DAT conversion
+    // Create a resampler for CD to DAT conversion (mono; for multi-channel
+    // audio use ProcessMulti, see "Multi-Channel Streaming" below)
     config := &resampling.Config{
         InputRate:  44100,
         OutputRate: 48000,
-        Channels:   2,
+        Channels:   1,
         Quality:    resampling.QualitySpec{Preset: resampling.QualityHigh},
     }
 
@@ -134,12 +135,49 @@ func main() {
         writeOutput(output)
     }
 
-    // Flush remaining samples
-    final, _ := r.Flush()
+    // Flush remaining samples at end of stream
+    final, err := r.Flush()
+    if err != nil {
+        log.Fatal(err)
+    }
     writeOutput(final)
 }
 ```
 
+### Latency and Real-Time Streaming
+
+A streaming resampler has a startup deficit: the internal filter needs a few samples of history before it can emit correctly filtered output, so the first `Process` calls in a stream withhold roughly `Latency()` samples that later calls make up. Callers that need a fixed number of output samples per callback, such as a portaudio or miniaudio audio callback, should sit a small FIFO between the resampler and the callback, primed with `Latency()` samples of silence. `Latency()` matches the measured deficit to within about 2 samples, so priming with it keeps callbacks fed in practice; any 1-2 sample shortfall self-heals because the deficit-driven buffer sizing catches up within the first few callbacks.
+
+`Latency()` is available on `SimpleResampler` and `SimpleResamplerFloat32` (the `NewEngine`/`NewEngineFloat32` path) only. Resamplers built from `New(config)` instead expose `GetLatency()`/`GetInfo()`, which report the filter group delay in the input domain: a different figure, not intended for FIFO priming.
+
+```go
+r, err := resampling.NewEngineFloat32(44100, 48000, resampling.QualityHigh)
+if err != nil {
+    log.Fatal(err)
+}
+fifo := make([]float32, r.Latency()) // prime with the startup deficit
+
+for chunk := range audioChunks {
+    out, err := r.Process(chunk)
+    if err != nil {
+        log.Fatal(err)
+    }
+    fifo = append(fifo, out...)
+    // deliver fixed-size slices from fifo to the audio callback here
+}
+
+// End of stream: drain the filter tail exactly once.
+tail, err := r.Flush()
+if err != nil {
+    log.Fatal(err)
+}
+fifo = append(fifo, tail...)
+```
+
+`Flush` is end-of-stream only. It pushes padding through the filter to drain its tail, so calling it once per chunk instead of once at the very end injects that padding into the middle of the stream and produces audible clicks at every chunk boundary. Call `Flush` exactly once, after the last `Process` call for the stream; a second `Flush` call is a no-op that returns an empty slice. `Reset` discards all filter state and is for starting an unrelated stream, never for use between chunks of the same stream. Each `SimpleResampler` or `SimpleResamplerFloat32` processes one channel, so a multi-channel stream keeps one persistent instance per channel alive for the duration of the stream.
+
+See [`examples/streaming`](examples/streaming/main.go) for the complete, runnable FIFO pattern.
+
 ### Zero-Allocation Streaming (`ProcessInto`)
 
 For allocation-sensitive pipelines, use caller-owned output buffers:
diff --git a/aliasing_test.go b/aliasing_test.go
new file mode 100644
index 0000000..cead8b7
--- /dev/null
+++ b/aliasing_test.go
@@ -0,0 +1,35 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import "testing"
+
+// Process must return an owned buffer at every ratio. At 1:1 the DFT stage
+// passthrough used to return the caller's own slice, so mutating the input
+// buffer afterwards corrupted previously returned output.
+func TestProcessOutputOwned_UnityRatio(t *testing.T) {
+	r, err := NewEngine(48000, 48000, QualityHigh)
+	if err != nil {
+		t.Fatal(err)
+	}
+	in := make([]float64, 256)
+	for i := range in {
+		in[i] = float64(i)
+	}
+	out, err := r.Process(in)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if len(out) != len(in) {
+		t.Fatalf("unity ratio length %d != %d", len(out), len(in))
+	}
+	for i := range in {
+		in[i] = -1
+	}
+	for i, v := range out {
+		if v != float64(i) {
+			t.Fatalf("output aliases input: out[%d] = %g after caller mutation", i, v)
+		}
+	}
+}
diff --git a/convenience.go b/convenience.go
index a4654d2..1ac2594 100644
--- a/convenience.go
+++ b/convenience.go
@@ -115,6 +115,10 @@ func NewMultiChannel(inputRate, outputRate float64, channels int, quality Qualit
 // SimpleResampler provides a simplified interface for basic resampling tasks.
 // It wraps the engine.Resampler directly for maximum performance.
 // Uses float64 precision for maximum quality.
+//
+// Streaming pattern: keep one persistent SimpleResampler per channel, call
+// Process once per chunk of input, and call Flush once at the end of the
+// stream.
 type SimpleResampler struct {
 	engine *engine.Resampler[float64]
 }
@@ -131,7 +135,12 @@ func NewEngine(inputRate, outputRate float64, quality QualityPreset) (*SimpleRes
 	return &SimpleResampler{engine: r}, nil
 }
 
-// Process resamples the input samples.
+// Process resamples one chunk of a single (mono) audio channel.
+// The output length varies from call to call: early calls may withhold up
+// to Latency() samples while the internal filter primes, so len(output)
+// does not track len(input)*GetRatio() on a per-call basis. The returned
+// slice is owned by the caller; it is never aliased by the resampler and
+// stays valid after subsequent calls.
 func (r *SimpleResampler) Process(input []float64) ([]float64, error) {
 	return r.engine.Process(input)
 }
@@ -165,17 +174,27 @@ func (r *SimpleResampler) EstimateOutput(inputLen int) int {
 	return int(float64(inputLen)*r.engine.GetRatio()) + estimateOutputMargin
 }
 
-// Flush returns any remaining buffered samples.
+// Flush drains the remaining buffered samples at end-of-stream. It is
+// end-of-stream only: pushing padding through the filter on every chunk
+// instead of just the last one creates audible edge transients.
+// Flush is terminal. After it returns, the instance behaves like a fresh
+// one for output purposes: a second Flush call returns an empty slice, and
+// a subsequent Process call produces output bit-identical to a fresh
+// instance. GetStatistics counters are cumulative and are not reset by
+// Flush.
 func (r *SimpleResampler) Flush() ([]float64, error) {
 	return r.engine.Flush()
 }
 
-// Reset clears internal state.
+// Reset discards all filter state, returning the resampler to its
+// just-constructed condition. It is for starting an unrelated stream, not
+// for use between chunks of one continuous stream: calling it mid-stream
+// destroys the filter's history and causes audible clicks.
 func (r *SimpleResampler) Reset() {
 	r.engine.Reset()
 }
 
-// GetRatio returns the resampling ratio.
+// GetRatio returns the resampling ratio (outputRate / inputRate).
 func (r *SimpleResampler) GetRatio() float64 {
 	return r.engine.GetRatio()
 }
@@ -185,10 +204,20 @@ func (r *SimpleResampler) GetStatistics() map[string]int64 {
 	return r.engine.GetStatistics()
 }
 
+// Latency returns the resampler's startup deficit in output samples: the
+// number of samples early Process calls withhold while the internal filter
+// primes. Real-time users feeding fixed-size output buffers should prime
+// their FIFO with this many samples of silence.
+func (r *SimpleResampler) Latency() int {
+	return r.engine.Latency()
+}
+
 // presetToEngineQuality converts a QualityPreset to engine.Quality.
 func presetToEngineQuality(preset QualityPreset) engine.Quality {
 	switch preset {
-	case QualityQuick, QualityLow:
+	case QualityQuick:
+		return engine.QualityQuick
+	case QualityLow:
 		return engine.QualityLow
 	case QualityMedium:
 		return engine.QualityMedium
@@ -200,7 +229,13 @@ func presetToEngineQuality(preset QualityPreset) engine.Quality {
 }
 
 // ResampleMono is a convenience function for one-shot mono resampling.
-// It creates a resampler, processes the input, flushes, and returns the result.
+// It creates a resampler, processes the input, flushes, and returns the
+// result. Use this when:
+//   - You have the entire input available upfront (not streaming)
+//   - float64 precision is required (mastering, archival)
+//
+// For real-time or chunked streaming, use NewEngine with Process per chunk
+// instead; see the streaming example.
 func ResampleMono(input []float64, inputRate, outputRate float64, quality QualityPreset) ([]float64, error) {
 	r, err := NewEngine(inputRate, outputRate, quality)
 	if err != nil {
@@ -301,6 +336,10 @@ func DeinterleaveFromStereo(interleaved []float64) (left, right []float64) {
 // returns float64), SimpleResamplerFloat32 keeps everything in float32,
 // eliminating type conversion overhead.
 //
+// Streaming pattern: keep one persistent SimpleResamplerFloat32 per channel,
+// call Process once per chunk of input, and call Flush once at the end of
+// the stream.
+//
 // Example:
 //
 //	r, err := resampler.NewEngineFloat32(44100, 48000, resampler.QualityHigh)
@@ -308,10 +347,17 @@ func DeinterleaveFromStereo(interleaved []float64) (left, right []float64) {
 //	    log.Fatal(err)
 //	}
 //	for chunk := range audioChunks {
-//	    output, _ := r.Process(chunk)  // []float32 in, []float32 out
+//	    output, err := r.Process(chunk)  // []float32 in, []float32 out
+//	    if err != nil {
+//	        log.Fatal(err)
+//	    }
 //	    writeOutput(output)
 //	}
-//	final, _ := r.Flush()  // Returns []float32!
+//	final, err := r.Flush()  // Returns []float32!
+//	if err != nil {
+//	    log.Fatal(err)
+//	}
+//	writeOutput(final)
 type SimpleResamplerFloat32 struct {
 	engine *engine.Resampler[float32]
 }
@@ -335,8 +381,13 @@ func NewEngineFloat32(inputRate, outputRate float64, quality QualityPreset) (*Si
 	return &SimpleResamplerFloat32{engine: r}, nil
 }
 
-// Process resamples the input samples.
+// Process resamples one chunk of a single (mono) audio channel.
 // Input and output are both float32, with no type conversion overhead.
+// The output length varies from call to call: early calls may withhold up
+// to Latency() samples while the internal filter primes, so len(output)
+// does not track len(input)*GetRatio() on a per-call basis. The returned
+// slice is owned by the caller; it is never aliased by the resampler and
+// stays valid after subsequent calls.
 func (r *SimpleResamplerFloat32) Process(input []float32) ([]float32, error) {
 	return r.engine.Process(input)
 }
@@ -372,14 +423,25 @@ func (r *SimpleResamplerFloat32) EstimateOutput(inputLen int) int {
 	return int(float64(inputLen)*r.engine.GetRatio()) + estimateOutputMargin
 }
 
-// Flush returns any remaining buffered samples as float32.
+// Flush drains the remaining buffered samples as float32 at end-of-stream.
 // Unlike the main Resampler.Flush() which returns float64, this returns
 // float32 for a consistent float32 workflow.
+//
+// Flush is end-of-stream only: pushing padding through the filter on every
+// chunk instead of just the last one creates audible edge transients (issue
+// #51). Flush is terminal. After it returns, the instance behaves like a
+// fresh one for output purposes: a second Flush call returns an empty
+// slice, and a subsequent Process call produces output bit-identical to a
+// fresh instance. GetStatistics counters are cumulative and are not reset
+// by Flush.
 func (r *SimpleResamplerFloat32) Flush() ([]float32, error) {
 	return r.engine.Flush()
 }
 
-// Reset clears internal state, allowing the resampler to be reused.
+// Reset discards all filter state, returning the resampler to its
+// just-constructed condition. It is for starting an unrelated stream, not
+// for use between chunks of one continuous stream: calling it mid-stream
+// destroys the filter's history and causes audible clicks.
 func (r *SimpleResamplerFloat32) Reset() {
 	r.engine.Reset()
 }
@@ -394,6 +456,14 @@ func (r *SimpleResamplerFloat32) GetStatistics() map[string]int64 {
 	return r.engine.GetStatistics()
 }
 
+// Latency returns the resampler's startup deficit in output samples: the
+// number of samples early Process calls withhold while the internal filter
+// primes. Real-time users feeding fixed-size output buffers should prime
+// their FIFO with this many samples of silence.
+func (r *SimpleResamplerFloat32) Latency() int {
+	return r.engine.Latency()
+}
+
 // ResampleMonoFloat32 is a convenience function for one-shot mono resampling
 // with float32 samples. It creates a resampler, processes the input, flushes,
 // and returns the result.
@@ -401,7 +471,10 @@ func (r *SimpleResamplerFloat32) GetStatistics() map[string]int64 {
 // This is the float32 equivalent of ResampleMono. Use this when:
 //   - Your audio data is already in float32 format
 //   - You want ~2x SIMD throughput compared to float64
-//   - 32-bit precision is sufficient (most real-time applications)
+//   - You have the entire input available upfront (not streaming)
+//
+// For real-time or chunked streaming, use NewEngineFloat32 with Process per
+// chunk instead; see the streaming example.
 //
 // For maximum precision (mastering, archival), use ResampleMono instead.
 func ResampleMonoFloat32(input []float32, inputRate, outputRate float64, quality QualityPreset) ([]float32, error) {
diff --git a/convenience_float32_test.go b/convenience_float32_test.go
index f16d255..2516b45 100644
--- a/convenience_float32_test.go
+++ b/convenience_float32_test.go
@@ -76,15 +76,18 @@ func TestSimpleResamplerFloat32_Process(t *testing.T) {
 	}
 }
 
-// TestSimpleResamplerFloat32_Flush verifies that Flush returns float32.
+// TestSimpleResamplerFloat32_Flush verifies that Flush returns the buffered
+// tail as finite float32 samples after a warm Process() call.
 func TestSimpleResamplerFloat32_Flush(t *testing.T) {
 	r, err := NewEngineFloat32(44100, 48000, QualityHigh)
 	if err != nil {
 		t.Fatalf("NewEngineFloat32 failed: %v", err)
 	}
 
-	// Process some samples first
-	input := make([]float32, 1000)
+	// 4410 samples (0.1s) is well beyond QualityHigh's filter latency, so
+	// the delay line is guaranteed to hold a real buffered tail; Flush must
+	// drain it rather than return empty.
+	input := make([]float32, 4410)
 	for i := range input {
 		input[i] = float32(math.Sin(2 * math.Pi * 1000 * float64(i) / 44100))
 	}
@@ -99,9 +102,13 @@ func TestSimpleResamplerFloat32_Flush(t *testing.T) {
 		t.Fatalf("Flush failed: %v", err)
 	}
 
-	// Flushed samples should exist (filter has latency)
 	if len(flushed) == 0 {
-		t.Log("Flush returned empty (may be valid depending on filter design)")
+		t.Fatal("Flush returned no samples after a warm Process() call; expected the buffered tail to drain")
+	}
+	for i, v := range flushed {
+		if math.IsNaN(float64(v)) || math.IsInf(float64(v), 0) {
+			t.Fatalf("flushed[%d] = %v, want finite", i, v)
+		}
 	}
 }
 
diff --git a/doc.go b/doc.go
index 2272bf4..d5f241d 100644
--- a/doc.go
+++ b/doc.go
@@ -37,12 +37,15 @@
 //	    log.Fatal(err)
 //	}
 //
-// For streaming resampling with a reusable resampler:
+// For streaming resampling with a reusable resampler (one mono channel; for
+// multi-channel audio call [Resampler.ProcessMulti] per chunk and, at
+// end-of-stream, [MultiFlusher.FlushMulti] via a type assertion on the
+// resampler, as shown in the [MultiFlusher] example):
 //
 //	config := &resampler.Config{
 //	    InputRate:  44100,
 //	    OutputRate: 48000,
-//	    Channels:   2,
+//	    Channels:   1,
 //	    Quality:    resampler.QualitySpec{Preset: resampler.QualityHigh},
 //	}
 //	r, err := resampler.New(config)
@@ -59,8 +62,12 @@
 //	    writeOutput(output)
 //	}
 //
-//	// Flush remaining samples
-//	final, _ := r.Flush()
+//	// Flush remaining samples at end of stream
+//	final, err := r.Flush()
+//	if err != nil {
+//	    log.Fatal(err)
+//	}
+//	writeOutput(final)
 //
 // # Zero-Allocation Streaming
 //
@@ -91,6 +98,48 @@
 // which run on the float32-native engine and are likewise zero-allocation once
 // warm.
 //
+// # Latency and Real-Time Streaming
+//
+// A streaming resampler has a startup deficit: the internal filter needs a
+// few samples of history before it can emit correctly filtered output, so
+// the first Process calls in a stream withhold roughly [SimpleResampler.Latency]
+// (or [SimpleResamplerFloat32.Latency]) samples that later calls make up.
+// Callers that need a fixed number of output samples per callback, such as
+// a portaudio or miniaudio audio callback, should sit a small FIFO between
+// the resampler and the callback, primed with Latency() samples of silence.
+// Latency() matches the measured deficit to within about 2 samples, so
+// priming with it keeps callbacks fed in practice; any 1-2 sample shortfall
+// self-heals as the deficit-driven buffer sizing catches up:
+//
+//	r, err := resampler.NewEngineFloat32(44100, 48000, resampler.QualityHigh)
+//	if err != nil {
+//	    log.Fatal(err)
+//	}
+//	fifo := make([]float32, r.Latency()) // prime with the startup deficit
+//	for chunk := range audioChunks {
+//	    out, err := r.Process(chunk)
+//	    if err != nil {
+//	        log.Fatal(err)
+//	    }
+//	    fifo = append(fifo, out...)
+//	    // deliver fixed-size slices from fifo to the audio callback here
+//	}
+//
+// [SimpleResampler.Flush] and [SimpleResamplerFloat32.Flush] are
+// end-of-stream only. Flush pushes padding through the filter to drain its
+// tail, so calling it once per chunk instead of once at the very end
+// injects that padding into the middle of the stream and produces audible
+// clicks at every chunk boundary. Call Flush exactly once, after the last
+// Process call for the stream; a second Flush call is a no-op that returns
+// an empty slice. [SimpleResampler.Reset] and [SimpleResamplerFloat32.Reset]
+// discard all filter state and are for starting an unrelated stream, never
+// for use between chunks of the same stream.
+//
+// Each SimpleResampler or SimpleResamplerFloat32 processes one channel, so
+// a multi-channel stream keeps one persistent instance per channel alive
+// for the duration of the stream. See examples/streaming for the complete,
+// runnable FIFO pattern.
+//
 // # Quality Presets
 //
 // The library provides several quality presets for common use cases:
@@ -137,10 +186,17 @@
 //	    log.Fatal(err)
 //	}
 //	for chunk := range audioChunks {
-//	    output, _ := r.Process(chunk)  // []float32 in, []float32 out
+//	    output, err := r.Process(chunk)  // []float32 in, []float32 out
+//	    if err != nil {
+//	        log.Fatal(err)
+//	    }
 //	    writeOutput(output)
 //	}
-//	final, _ := r.Flush()  // Returns []float32 (not []float64!)
+//	final, err := r.Flush()  // Returns []float32 (not []float64!)
+//	if err != nil {
+//	    log.Fatal(err)
+//	}
+//	writeOutput(final)
 //
 // Helper functions for float32 stereo interleaving are also provided:
 // [InterleaveToStereoFloat32] and [DeinterleaveFromStereoFloat32].
diff --git a/examples/streaming/main.go b/examples/streaming/main.go
new file mode 100644
index 0000000..3a921db
--- /dev/null
+++ b/examples/streaming/main.go
@@ -0,0 +1,110 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+// Example: real-time chunked resampling into fixed-size output buffers.
+//
+// An audio callback (portaudio, miniaudio, ...) demands exactly N output
+// frames per call, but a streaming resampler returns a varying number of
+// samples per Process call (early calls withhold samples while the filter
+// primes). The fix is a small FIFO between the resampler and the callback,
+// primed with Latency() samples of silence. Never call Flush or Reset
+// inside the stream: Flush is end-of-stream only, Reset starts a new
+// stream (both destroy continuity and cause audible clicks, issue #51).
+package main
+
+import (
+	"fmt"
+	"math"
+
+	resampler "github.com/tphakala/go-audio-resampler"
+)
+
+func main() {
+	const (
+		inRate        = 44100.0
+		outRate       = 48000.0
+		outFrames     = 512
+		callbacks     = 100
+		toneAmplitude = 0.5
+		toneHz        = 997.0
+	)
+
+	rs, err := resampler.NewEngineFloat32(inRate, outRate, resampler.QualityHigh)
+	if err != nil {
+		panic(err)
+	}
+	ratio := rs.GetRatio()
+
+	// Prime the FIFO with the startup deficit so the first callbacks are
+	// fed. This trades Latency() samples of leading silence for a steady
+	// pipeline. The FIFO starts empty with headroom, then the priming zeros
+	// are appended so later appends grow a zero-length-origin slice.
+	fifo := make([]float32, 0, rs.Latency()+2*outFrames)
+	fifo = append(fifo, make([]float32, rs.Latency())...)
+
+	phase := 0.0
+	firstCall := true
+	for callback := range callbacks {
+		// Size the input chunk from the FIFO's current deficit rather than
+		// a fixed count. A fixed input size drifts against a fixed output
+		// size whenever ratio does not divide outFrames evenly: truncating
+		// the fixed size underfeeds the resampler and causes underruns
+		// every few seconds, while rounding it up overfeeds and grows the
+		// FIFO (and its latency) without bound over a long-running stream.
+		// Pulling exactly enough input to cover the current shortfall self-
+		// corrects both directions and keeps the FIFO bounded.
+		//
+		// The very first Process call is a special case: a fresh engine
+		// pays its entire Latency() startup deficit on that one call,
+		// regardless of how much input it receives, and the priming above
+		// exists to cover exactly that. So the first request must target a
+		// full outFrames, not outFrames minus the priming already sitting
+		// in the FIFO; netting the priming against the first request would
+		// count the same deficit twice and under-deliver on that call.
+		//
+		// This holds when outFrames comfortably exceeds Latency(), as here.
+		// With an output buffer smaller than the deficit, a single primed
+		// Process cannot cover it and the shortfall instead spreads over the
+		// first several callbacks until the FIFO fills; the underrun branch
+		// below tolerates that warmup.
+		need := outFrames
+		if !firstCall {
+			need = max(outFrames-len(fifo), 0)
+		}
+		firstCall = false
+		inFrames := int(math.Ceil(float64(need) / ratio))
+
+		// Allocated per callback for clarity; a production callback should
+		// reuse a single scratch buffer instead of allocating each call.
+		in := make([]float32, inFrames)
+		for i := range in {
+			in[i] = float32(toneAmplitude * math.Sin(phase))
+			phase += 2 * math.Pi * toneHz / inRate
+		}
+		out, err := rs.Process(in)
+		if err != nil {
+			panic(err)
+		}
+		fifo = append(fifo, out...)
+
+		if len(fifo) >= outFrames {
+			deliver := fifo[:outFrames]
+			_ = deliver // hand exactly outFrames samples to the audio API here
+			// Resliced from the front for clarity; the drained head keeps the
+			// backing array growing over a long stream, so production code
+			// would use a ring buffer instead.
+			fifo = fifo[outFrames:]
+		} else {
+			// Underrun (should not happen after priming): deliver silence.
+			fmt.Printf("callback %d: FIFO underrun (%d < %d)\n", callback, len(fifo), outFrames)
+		}
+	}
+
+	// End of stream: drain the filter tail exactly once.
+	tail, err := rs.Flush()
+	if err != nil {
+		panic(err)
+	}
+	fifo = append(fifo, tail...)
+	fmt.Printf("stream done, %d samples left to deliver\n", len(fifo))
+}
diff --git a/flush_length_test.go b/flush_length_test.go
new file mode 100644
index 0000000..5c50fe1
--- /dev/null
+++ b/flush_length_test.go
@@ -0,0 +1,54 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import (
+	"math"
+	"testing"
+)
+
+// Process+Flush of a fresh instance must emit exactly the resampled length:
+// no phantom padding samples. Issue #51: 470 samples at 44100 to 48000
+// QualityHigh returned 514 where ceil(470*ratio)+1 = 513 is the maximum.
+func TestFlushLength_Canonical(t *testing.T) {
+	cases := []struct {
+		in, out float64
+		n       int
+	}{
+		{44100, 48000, 470},
+		{44100, 48000, 4410},
+		{44100, 48000, 44100},
+		{48000, 44100, 480},
+		{48000, 44100, 48000},
+		{48000, 16000, 4800},
+	}
+	for _, q := range []QualityPreset{QualityMedium, QualityHigh} {
+		for _, c := range cases {
+			r, err := NewEngine(c.in, c.out, q)
+			if err != nil {
+				t.Fatal(err)
+			}
+			in := make([]float64, c.n)
+			for i := range in {
+				in[i] = 0.5 * math.Sin(2*math.Pi*997*float64(i)/c.in)
+			}
+			out, err := r.Process(in)
+			if err != nil {
+				t.Fatal(err)
+			}
+			tail, err := r.Flush()
+			if err != nil {
+				t.Fatal(err)
+			}
+			total := len(out) + len(tail)
+			ideal := float64(c.n) * c.out / c.in
+			lo := int(math.Floor(ideal))
+			hi := int(math.Ceil(ideal)) + 1
+			if total < lo || total > hi {
+				t.Errorf("q=%v %v to %v n=%d: total %d outside [%d, %d]",
+					q, c.in, c.out, c.n, total, lo, hi)
+			}
+		}
+	}
+}
diff --git a/flush_lifecycle_test.go b/flush_lifecycle_test.go
new file mode 100644
index 0000000..806e712
--- /dev/null
+++ b/flush_lifecycle_test.go
@@ -0,0 +1,78 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import (
+	"math"
+	"slices"
+	"testing"
+)
+
+func sineChunk(n int, rate float64) []float64 {
+	out := make([]float64, n)
+	for i := range out {
+		out[i] = 0.5 * math.Sin(2*math.Pi*997*float64(i)/rate)
+	}
+	return out
+}
+
+// Flush ends the stream: a second Flush must return nothing, and a
+// subsequent Process must behave exactly like a fresh instance instead of
+// convolving new audio against leftover padding zeros.
+func TestFlushLifecycle(t *testing.T) {
+	for _, c := range []struct {
+		in, out float64
+		quality QualityPreset
+	}{
+		{44100, 48000, QualityHigh},
+		{48000, 44100, QualityHigh},
+		{48000, 16000, QualityHigh},
+		// QualityQuick routes through CubicStage (reachable from NewEngine
+		// since the QualityQuick mapping fix), which has its own held-tail
+		// lifecycle distinct from the FIR stages' delay lines.
+		{44100, 48000, QualityQuick},
+	} {
+		r, err := NewEngine(c.in, c.out, c.quality)
+		if err != nil {
+			t.Fatal(err)
+		}
+		if _, err := r.Process(sineChunk(4410, c.in)); err != nil {
+			t.Fatal(err)
+		}
+		if _, err := r.Flush(); err != nil {
+			t.Fatal(err)
+		}
+
+		second, err := r.Flush()
+		if err != nil {
+			t.Fatal(err)
+		}
+		if len(second) != 0 {
+			t.Errorf("%v to %v q=%v: second Flush returned %d samples, want 0", c.in, c.out, c.quality, len(second))
+		}
+
+		fresh, err := NewEngine(c.in, c.out, c.quality)
+		if err != nil {
+			t.Fatal(err)
+		}
+		chunk := sineChunk(4410, c.in)
+		gotAfterFlush, err := r.Process(slices.Clone(chunk))
+		if err != nil {
+			t.Fatal(err)
+		}
+		gotFresh, err := fresh.Process(slices.Clone(chunk))
+		if err != nil {
+			t.Fatal(err)
+		}
+		if len(gotAfterFlush) != len(gotFresh) {
+			t.Fatalf("%v to %v q=%v: post-flush Process length %d != fresh %d",
+				c.in, c.out, c.quality, len(gotAfterFlush), len(gotFresh))
+		}
+		for i := range gotFresh {
+			if gotAfterFlush[i] != gotFresh[i] {
+				t.Fatalf("%v to %v q=%v: post-flush Process differs from fresh at %d", c.in, c.out, c.quality, i)
+			}
+		}
+	}
+}
diff --git a/flush_multi_test.go b/flush_multi_test.go
index 27c3482..031ad7d 100644
--- a/flush_multi_test.go
+++ b/flush_multi_test.go
@@ -7,6 +7,7 @@ import (
 	"fmt"
 	"math"
 	"math/rand"
+	"slices"
 	"testing"
 )
 
@@ -98,8 +99,11 @@ func TestFlushMulti_MatchesPerChannelFlush(t *testing.T) {
 		channelInputs[ch] = makeDeterministicInput(rng, n, inRate)
 	}
 
-	// Per-channel reference: process each channel independently with a mono resampler.
-	monoTotals := make([]int, channels)
+	// Per-channel reference: process each channel independently with a mono
+	// resampler, keeping the full Process+Flush sample content (not just its
+	// length) so the multi-channel path can be checked for content, not
+	// merely count.
+	monoOutputs := make([][]float64, channels)
 	for ch := range channels {
 		cfg := &Config{
 			InputRate:  inRate,
@@ -119,7 +123,7 @@ func TestFlushMulti_MatchesPerChannelFlush(t *testing.T) {
 		if err != nil {
 			t.Fatalf("mono Flush ch%d: %v", ch, err)
 		}
-		monoTotals[ch] = len(proc) + len(fl)
+		monoOutputs[ch] = append(slices.Clone(proc), fl...)
 	}
 
 	// Multi-channel: process all channels together.
@@ -148,10 +152,19 @@ func TestFlushMulti_MatchesPerChannelFlush(t *testing.T) {
 	}
 
 	for ch := range channels {
-		multiTotal := len(proc[ch]) + len(flushed[ch])
-		if multiTotal != monoTotals[ch] {
-			t.Errorf("channel %d: multi total %d != mono total %d",
-				ch, multiTotal, monoTotals[ch])
+		multiOutput := append(slices.Clone(proc[ch]), flushed[ch]...)
+		if len(multiOutput) != len(monoOutputs[ch]) {
+			t.Fatalf("channel %d: multi total %d != mono total %d",
+				ch, len(multiOutput), len(monoOutputs[ch]))
+		}
+		// Fatalf (not Errorf) stops at the first mismatch: a real regression
+		// here differs across most of the signal, and letting the loop run
+		// to completion would flood the test log with thousands of lines.
+		for i := range multiOutput {
+			if multiOutput[i] != monoOutputs[ch][i] {
+				t.Fatalf("channel %d: sample %d differs: multi=%v mono=%v",
+					ch, i, multiOutput[i], monoOutputs[ch][i])
+			}
 		}
 	}
 }
@@ -180,4 +193,9 @@ func TestFlushMulti_EmptyResampler(t *testing.T) {
 	if len(flushed) != 2 {
 		t.Fatalf("FlushMulti returned %d channels, want 2", len(flushed))
 	}
+	for ch, samples := range flushed {
+		if len(samples) != 0 {
+			t.Errorf("channel %d: FlushMulti on an empty resampler returned %d samples, want 0", ch, len(samples))
+		}
+	}
 }
diff --git a/internal/engine/buffer_integrity_test.go b/internal/engine/buffer_integrity_test.go
index 282e947..03613e1 100644
--- a/internal/engine/buffer_integrity_test.go
+++ b/internal/engine/buffer_integrity_test.go
@@ -283,82 +283,125 @@ func TestResampler_ProcessAndFlushSequence(t *testing.T) {
 }
 
 // TestDFTStage_MultipleProcessCalls tests that multiple consecutive Process()
-// calls work correctly and don't accumulate errors.
+// calls work correctly, don't corrupt earlier returned outputs, and are
+// deterministic: an identical call sequence on a second fresh instance must
+// reproduce the same outputs bit-exactly.
 func TestDFTStage_MultipleProcessCalls(t *testing.T) {
-	stage, err := NewDFTStage[float64](2, QualityHigh)
-	require.NoError(t, err, "Failed to create DFT stage")
-
 	const numCalls = 10
 	const samplesPerCall = 500
 
-	// Store all outputs
-	outputs := make([][]float64, numCalls)
-
-	for i := range numCalls {
+	// makeInput builds the same deterministic per-call input for both
+	// instances below, so the two runs are directly comparable.
+	makeInput := func(callIdx int) []float64 {
 		input := make([]float64, samplesPerCall)
 		for j := range input {
 			// Different phase for each call to detect any cross-contamination
-			input[j] = math.Sin(2.0*math.Pi*float64(j)/100 + float64(i)*math.Pi/5)
+			input[j] = math.Sin(2.0*math.Pi*float64(j)/100 + float64(callIdx)*math.Pi/5)
+		}
+		return input
+	}
+
+	runSequence := func(stage *DFTStage[float64]) (raw, saved [][]float64) {
+		raw = make([][]float64, numCalls)
+		saved = make([][]float64, numCalls)
+		for i := range numCalls {
+			output, err := stage.Process(makeInput(i))
+			require.NoError(t, err, "Process() call %d failed", i)
+			raw[i] = output
+			saved[i] = make([]float64, len(output))
+			copy(saved[i], output)
 		}
+		return raw, saved
+	}
 
-		output, err := stage.Process(input)
-		require.NoError(t, err, "Process() call %d failed", i)
+	stageA, err := NewDFTStage[float64](2, QualityHigh)
+	require.NoError(t, err, "Failed to create DFT stage")
+	rawA, savedA := runSequence(stageA)
 
-		// Save a copy
-		outputs[i] = make([]float64, len(output))
-		copy(outputs[i], output)
+	// Every earlier call's returned slice must still hold the values it held
+	// right after that call: a later Process() reusing the same backing
+	// array without properly copying out would silently corrupt it.
+	for i := range numCalls {
+		require.Len(t, rawA[i], len(savedA[i]), "call %d output length changed", i)
+		for j := range rawA[i] {
+			// require (not assert) stops at the first mismatch: a real
+			// corruption bug can differ across most of a call's samples,
+			// and letting the loop run to completion floods the test log.
+			require.InDelta(t, savedA[i][j], rawA[i][j], 1e-15,
+				"call %d output[%d] was corrupted by a later Process() call", i, j)
+		}
 	}
 
-	// Verify all stored outputs are still valid (not corrupted by subsequent calls)
-	for i := range numCalls - 1 {
-		for j := range min(10, len(outputs[i])) {
-			// The saved value should match what we stored
-			assert.False(t, math.IsNaN(outputs[i][j]),
-				"outputs[%d][%d] became NaN", i, j)
-			assert.False(t, math.IsInf(outputs[i][j], 0),
-				"outputs[%d][%d] became Inf", i, j)
+	// Determinism: an identical call sequence on a second, fresh instance
+	// must produce bit-identical output per call.
+	stageB, err := NewDFTStage[float64](2, QualityHigh)
+	require.NoError(t, err, "Failed to create second DFT stage")
+	_, savedB := runSequence(stageB)
+
+	for i := range numCalls {
+		require.Len(t, savedB[i], len(savedA[i]), "call %d length differs between two fresh instances", i)
+		for j := range savedA[i] {
+			require.InDelta(t, savedA[i][j], savedB[i][j], 1e-15,
+				"call %d output[%d] differs between two fresh instances given identical input", i, j)
 		}
 	}
 
-	t.Logf("DFT stage: %d consecutive Process() calls verified", numCalls)
+	t.Logf("DFT stage: %d consecutive Process() calls verified bit-identical and uncorrupted", numCalls)
 }
 
-// TestPolyphaseStage_MultipleProcessCalls tests that multiple consecutive Process()
-// calls work correctly.
+// TestPolyphaseStage_MultipleProcessCalls tests that multiple consecutive
+// Process() calls work correctly, don't corrupt earlier returned outputs,
+// and are deterministic across a second fresh instance.
 func TestPolyphaseStage_MultipleProcessCalls(t *testing.T) {
-	stage, err := NewPolyphaseStage[float64](1.088435374, 0.459375, true, QualityHigh)
-	require.NoError(t, err, "Failed to create polyphase stage")
-
 	const numCalls = 10
 	const samplesPerCall = 1000
 
-	// Store all outputs
-	outputs := make([][]float64, numCalls)
-
-	for i := range numCalls {
+	makeInput := func(callIdx int) []float64 {
 		input := make([]float64, samplesPerCall)
 		for j := range input {
-			input[j] = math.Sin(2.0*math.Pi*float64(j)/100 + float64(i)*math.Pi/5)
+			input[j] = math.Sin(2.0*math.Pi*float64(j)/100 + float64(callIdx)*math.Pi/5)
+		}
+		return input
+	}
+
+	runSequence := func(stage *PolyphaseStage[float64]) (raw, saved [][]float64) {
+		raw = make([][]float64, numCalls)
+		saved = make([][]float64, numCalls)
+		for i := range numCalls {
+			output, err := stage.Process(makeInput(i))
+			require.NoError(t, err, "Process() call %d failed", i)
+			raw[i] = output
+			saved[i] = make([]float64, len(output))
+			copy(saved[i], output)
 		}
+		return raw, saved
+	}
 
-		output, err := stage.Process(input)
-		require.NoError(t, err, "Process() call %d failed", i)
+	stageA, err := NewPolyphaseStage[float64](1.088435374, 0.459375, true, QualityHigh)
+	require.NoError(t, err, "Failed to create polyphase stage")
+	rawA, savedA := runSequence(stageA)
 
-		outputs[i] = make([]float64, len(output))
-		copy(outputs[i], output)
+	for i := range numCalls {
+		require.Len(t, rawA[i], len(savedA[i]), "call %d output length changed", i)
+		for j := range rawA[i] {
+			require.InDelta(t, savedA[i][j], rawA[i][j], 1e-15,
+				"call %d output[%d] was corrupted by a later Process() call", i, j)
+		}
 	}
 
-	// Verify all stored outputs are still valid
-	for i := range numCalls - 1 {
-		for j := range min(10, len(outputs[i])) {
-			assert.False(t, math.IsNaN(outputs[i][j]),
-				"outputs[%d][%d] became NaN", i, j)
-			assert.False(t, math.IsInf(outputs[i][j], 0),
-				"outputs[%d][%d] became Inf", i, j)
+	stageB, err := NewPolyphaseStage[float64](1.088435374, 0.459375, true, QualityHigh)
+	require.NoError(t, err, "Failed to create second polyphase stage")
+	_, savedB := runSequence(stageB)
+
+	for i := range numCalls {
+		require.Len(t, savedB[i], len(savedA[i]), "call %d length differs between two fresh instances", i)
+		for j := range savedA[i] {
+			require.InDelta(t, savedA[i][j], savedB[i][j], 1e-15,
+				"call %d output[%d] differs between two fresh instances given identical input", i, j)
 		}
 	}
 
-	t.Logf("Polyphase stage: %d consecutive Process() calls verified", numCalls)
+	t.Logf("Polyphase stage: %d consecutive Process() calls verified bit-identical and uncorrupted", numCalls)
 }
 
 // TestCubicStage_BufferIntegrity verifies CubicStage doesn't have buffer issues.
diff --git a/internal/engine/cubic.go b/internal/engine/cubic.go
index 03447e3..ccc2830 100644
--- a/internal/engine/cubic.go
+++ b/internal/engine/cubic.go
@@ -13,11 +13,12 @@ import (
 // CubicStage implements cubic (4-point, 3rd order) interpolation matching SOXR.
 // This is the fastest resampling method, used for QualityQuick preset.
 type CubicStage[F simdops.Float] struct {
-	ratio   float64
-	phase   float64
-	history [4]F // 4-point window for interpolation
-	histPos int
-	latency int
+	ratio     float64
+	phase     float64
+	history   [4]F // 4-point window for interpolation
+	primed    int  // real samples pushed so far, capped at cubicLatencySamples
+	latency   int
+	outputBuf []F // reused across calls so the warm path allocates nothing
 }
 
 // NewCubicStage creates a new cubic interpolation stage.
@@ -29,15 +30,46 @@ func NewCubicStage[F simdops.Float](ratio float64) *CubicStage[F] {
 	}
 }
 
-// Process resamples input using cubic interpolation.
+// Process resamples input using cubic interpolation. The returned slice is
+// owned by the caller and remains valid across subsequent calls.
 func (c *CubicStage[F]) Process(input []F) ([]F, error) {
+	out, err := c.processZeroCopy(input)
+	if err != nil {
+		return out, err
+	}
+	if len(out) == 0 {
+		// processZeroCopy returns an empty slice that may still alias
+		// c.outputBuf (len 0, cap > 0) during the priming window. Return a
+		// fresh literal so a caller appending to the result is not corrupted
+		// when the next call reuses the internal buffer.
+		return []F{}, nil
+	}
+	// Return a copy so the caller's slice is not corrupted when the next call
+	// reuses the internal output buffer.
+	result := make([]F, len(out))
+	copy(result, out)
+	return result, nil
+}
+
+// processZeroCopy is the allocation-free internal path. The returned slice
+// aliases c.outputBuf and is only valid until the next Process,
+// processZeroCopy, Flush, or Reset call.
+func (c *CubicStage[F]) processZeroCopy(input []F) ([]F, error) { //nolint:unparam // error kept for symmetry with the FIR stages' Process signature
 	if len(input) == 0 {
 		return []F{}, nil
 	}
 
-	// Estimate output size
-	outputSize := int(math.Ceil(float64(len(input)) * c.ratio))
-	output := make([]F, 0, outputSize)
+	// Upper bound on the outputs this call can emit: the phase accumulator
+	// advances one input unit per sample and emits at most one output per
+	// 1/ratio of that advance, so the count never exceeds ceil(len*ratio)
+	// plus one boundary-carry sample. Pre-size the reused buffer to that
+	// bound via growStableLen (which adds headroom when it must grow) so
+	// steady-state constant-chunk streaming reuses it without allocating.
+	// append then fills it, so an off-by-one in the bound can never write out
+	// of range; on the warm path the capacity already covers the count and
+	// append allocates nothing.
+	maxOut := int(math.Ceil(float64(len(input))*c.ratio)) + 1
+	out := growStableLen(c.outputBuf, maxOut)[:0]
 
 	for _, sample := range input {
 		// Shift history window
@@ -46,11 +78,27 @@ func (c *CubicStage[F]) Process(input []F) ([]F, error) {
 		c.history[1] = c.history[0]
 		c.history[0] = sample
 
+		// The center point used below (history[2]) needs cubicLatencySamples
+		// real pushes past it before it holds real signal instead of the
+		// zero value left by construction or Reset. Withhold emission during
+		// that priming window: emission starts only once history[2] holds a
+		// real sample. The oldest neighbor (history[3]) can still be the zero
+		// initial state for the first post-priming emission, a bounded startup
+		// characteristic, not fabricated steady-state output. Flush drains the
+		// true tail this reserves.
+		//
+		// Invariant: cubicLatencySamples must equal the shift depth from the
+		// push at history[0] to the center at history[2] (two shifts); the
+		// priming gate relies on that coincidence.
+		if c.primed < cubicLatencySamples {
+			c.primed++
+			continue
+		}
+
 		// Generate output samples
 		for c.phase < 1.0 {
 			// Cubic interpolation matching SOXR
-			y := c.interpolate(c.phase)
-			output = append(output, y)
+			out = append(out, c.interpolate(c.phase))
 
 			// Advance phase
 			c.phase += 1.0 / c.ratio
@@ -60,7 +108,8 @@ func (c *CubicStage[F]) Process(input []F) ([]F, error) {
 		c.phase -= 1.0
 	}
 
-	return output, nil
+	c.outputBuf = out
+	return out, nil
 }
 
 // interpolate performs cubic interpolation matching SOXR's implementation.
@@ -89,16 +138,38 @@ func (c *CubicStage[F]) interpolate(x float64) F {
 	return F(((a*x+b)*x+coefC)*x + s0)
 }
 
-// Flush returns any remaining samples.
+// Flush drains the interpolator's held tail.
+//
+// The 4-point history window holds cubicLatencySamples of latency: the
+// final real samples pushed into Process never reach the withheld-emission
+// center position (see Process) and stay trapped internally when the caller
+// stops feeding input. Flush pads that many zeros through the same Process
+// path to release them, then resets to fresh state so a second Flush
+// returns nothing and a subsequent Process behaves like a new instance.
+//
+// The zero padding is a hard silence assumption at the true end of signal.
+// For a signal with a large sample-to-sample swing right at the boundary,
+// cubic's polynomial fit can overshoot past the padding transition before
+// settling; unlike an FIR filter's convolution, which decays smoothly by
+// construction, point interpolation has no equivalent damping. This is a
+// known, minor characteristic, not something Flush can avoid while still
+// delivering the real tail.
 func (c *CubicStage[F]) Flush() ([]F, error) {
-	// Cubic interpolation doesn't buffer samples
-	return []F{}, nil
+	if c.primed == 0 {
+		// Never fed: no held tail to drain.
+		return []F{}, nil
+	}
+	zeros := make([]F, cubicLatencySamples)
+	out, err := c.Process(zeros)
+	c.Reset()
+	return out, err
 }
 
 // Reset clears internal state.
 func (c *CubicStage[F]) Reset() {
 	c.phase = 0
 	c.history = [4]F{}
+	c.primed = 0
 }
 
 // GetRatio returns the stage's resampling ratio.
diff --git a/internal/engine/cubic_flush_test.go b/internal/engine/cubic_flush_test.go
new file mode 100644
index 0000000..a44202b
--- /dev/null
+++ b/internal/engine/cubic_flush_test.go
@@ -0,0 +1,334 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package engine
+
+import (
+	"math"
+	"slices"
+	"testing"
+)
+
+// The cubic interpolator holds a 4-point window with 2 samples of latency;
+// Flush must emit the tail those samples cover instead of dropping it.
+//
+// Total length alone is not a reliable regression signal for this stage:
+// CubicStage.Process (unlike the FIR stages) never withholds output pending
+// future context, so it emits immediately from a fictional zero-filled
+// history. That keeps Process-alone output count close to n*ratio even when
+// the true final samples are dropped, because emitting a few samples from
+// fictional pre-silence at the head happens to offset the count. So this
+// test also checks head content: the first few output samples must track
+// the ramp continuously instead of dipping toward zero from that fictional
+// history.
+func TestCubicStage_FlushEmitsTail(t *testing.T) {
+	r, err := NewResampler[float64](44100, 48000, QualityQuick)
+	if err != nil {
+		t.Fatal(err)
+	}
+	const n = 4410
+	in := make([]float64, n)
+	for i := range in {
+		in[i] = float64(i)
+	}
+	out, err := r.Process(in)
+	if err != nil {
+		t.Fatal(err)
+	}
+	tail, err := r.Flush()
+	if err != nil {
+		t.Fatal(err)
+	}
+	total := len(out) + len(tail)
+	ideal := int(float64(n) * 48000.0 / 44100.0)
+	if total < ideal-1 || total > ideal+1 {
+		t.Errorf("cubic Process+Flush total %d, want %d +-1", total, ideal)
+	}
+
+	if len(out) < 4 {
+		t.Fatalf("cubic Process produced only %d samples, too few to check head content", len(out))
+	}
+	const step = 44100.0 / 48000.0
+	for i := range 4 {
+		want := float64(i) * step
+		got := out[i]
+		if math.Abs(got-want) > 0.5 {
+			t.Errorf("head sample %d = %v, want within 0.5 of %v (ramp continuity, not fictional zero history)", i, got, want)
+		}
+	}
+}
+
+// The flushed tail must contain the real trailing segments (centered on
+// x[n-2] and x[n-1]), not garbage: its first sample should still track the
+// ramp's final value, with later samples decaying toward the zero padding.
+//
+// This deliberately does not reuse TestCubicStage_FlushEmitsTail's unbounded
+// 0..4409 ramp. That ramp is well suited to the head check (early values are
+// widely spaced integers, so a fictional-zero head is obviously distinct
+// from a correct one) but poorly suited to a tail check: its final value is
+// thousands of units away from the zero Flush pads in, and cubic
+// interpolation's polynomial fit overshoots substantially in the interior of
+// a segment spanning that large a discontinuity (its endpoints stay exact;
+// only the interior curve swings, by hundreds of units for that ramp's
+// scale) even on the already-fixed implementation. This is the same known,
+// minor characteristic documented on CubicStage.Flush, not a defect; a
+// bounded ramp keeps the discontinuity small enough that tail values are
+// checkable with a tight, meaningful tolerance instead of one loose enough
+// to hide a real regression.
+func TestCubicStage_FlushTailTracksRamp(t *testing.T) {
+	c := NewCubicStage[float64](48000.0 / 44100.0)
+	const n = 4410
+	in := make([]float64, n)
+	for i := range in {
+		in[i] = float64(i) / float64(n-1) // bounded 0..1 ramp
+	}
+	if _, err := c.Process(in); err != nil {
+		t.Fatal(err)
+	}
+	tail, err := c.Flush()
+	if err != nil {
+		t.Fatal(err)
+	}
+	if len(tail) == 0 {
+		t.Fatal("Flush returned no samples; nothing to check")
+	}
+
+	const tolerance = 0.1
+	last := in[n-1]
+	if math.Abs(tail[0]-last) > tolerance {
+		t.Errorf("tail[0] = %v, want within %v of ramp end %v (real tail, not dropped)", tail[0], tolerance, last)
+	}
+	if math.Abs(tail[len(tail)-1]) > tolerance {
+		t.Errorf("final tail sample = %v, want within %v of 0 (decayed toward the zero padding)", tail[len(tail)-1], tolerance)
+	}
+}
+
+// Process exactly one real sample (0 < primed < cubicLatencySamples), then
+// Flush. This is the partially-primed edge: the interpolator never reaches
+// the priming threshold during Process, so Flush's own zero padding must
+// both finish priming and drain the single real sample, without a panic and
+// without breaking the terminal-flush lifecycle.
+func TestCubicStage_FlushAfterPartialPriming(t *testing.T) {
+	c := NewCubicStage[float64](48000.0 / 44100.0)
+
+	out, err := c.Process([]float64{42.0})
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	tail, err := c.Flush()
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	total := len(out) + len(tail)
+	n := 1 // runtime variable: 1*48000/44100 isn't an exact integer constant
+	ideal := int(float64(n) * 48000.0 / 44100.0)
+	if total < ideal-1 || total > ideal+1 {
+		t.Errorf("partially primed Process+Flush total %d, want %d +-1", total, ideal)
+	}
+
+	// ideal is 1 here, so the +-1 count check above alone tolerates
+	// total=0: a regression that silently drops the single pushed sample
+	// would still pass it. Pin the drain directly: the tail must be
+	// non-empty, and cubic interpolation evaluates exactly to the center
+	// history point at x=0 (all polynomial terms multiply by x and vanish),
+	// so tail[0] must equal the pushed sample.
+	if len(tail) == 0 {
+		t.Fatal("Flush returned no samples for a single partially-primed input; the pushed sample was dropped")
+	}
+	const pushedSample = 42.0
+	const tolerance = 1e-9
+	if math.Abs(tail[0]-pushedSample) > tolerance {
+		t.Errorf("tail[0] = %v, want within %v of pushed sample %v", tail[0], tolerance, pushedSample)
+	}
+
+	second, err := c.Flush()
+	if err != nil {
+		t.Fatal(err)
+	}
+	if len(second) != 0 {
+		t.Errorf("second Flush returned %d samples, want 0", len(second))
+	}
+
+	fresh := NewCubicStage[float64](48000.0 / 44100.0)
+	freshOut, err := fresh.Process([]float64{42.0})
+	if err != nil {
+		t.Fatal(err)
+	}
+	gotAfterFlush, err := c.Process([]float64{42.0})
+	if err != nil {
+		t.Fatal(err)
+	}
+	if len(gotAfterFlush) != len(freshOut) {
+		t.Fatalf("post-flush Process length %d != fresh %d", len(gotAfterFlush), len(freshOut))
+	}
+	for i := range freshOut {
+		if gotAfterFlush[i] != freshOut[i] {
+			t.Fatalf("post-flush Process differs from fresh at %d", i)
+		}
+	}
+}
+
+// Flush must drain the true tail, not just leave the stage silently short:
+// after real input has primed the interpolator, Flush must return samples.
+func TestCubicStage_FlushIsNonEmptyAfterRealInput(t *testing.T) {
+	c := NewCubicStage[float64](48000.0 / 44100.0)
+	in := make([]float64, 100)
+	for i := range in {
+		in[i] = float64(i)
+	}
+	if _, err := c.Process(in); err != nil {
+		t.Fatal(err)
+	}
+	tail, err := c.Flush()
+	if err != nil {
+		t.Fatal(err)
+	}
+	if len(tail) == 0 {
+		t.Fatal("Flush returned no samples after real input was processed; held tail was dropped")
+	}
+}
+
+// GetStatistics()["samplesOut"] must count every sample that Process and
+// Flush return. On the FIR path the engine's Flush adds len(output) to
+// samplesOut; the cubic (QualityQuick) branch early-returns the stage's
+// flush tail and used to skip that accounting, so the statistic undercounted
+// by the flush-tail length once cubic Flush began emitting a real tail. This
+// pins the invariant: samplesOut == len(Process output) + len(Flush output).
+func TestCubicStage_FlushUpdatesSamplesOut(t *testing.T) {
+	r, err := NewResampler[float64](44100, 48000, QualityQuick)
+	if err != nil {
+		t.Fatal(err)
+	}
+	const n = 4410
+	in := make([]float64, n)
+	for i := range in {
+		in[i] = float64(i)
+	}
+	out, err := r.Process(in)
+	if err != nil {
+		t.Fatal(err)
+	}
+	tail, err := r.Flush()
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	// A non-empty tail is what makes the undercount observable; without it
+	// the assertion below could pass even with the accounting bug present.
+	if len(tail) == 0 {
+		t.Fatal("Flush returned no tail; test cannot distinguish the samplesOut undercount")
+	}
+
+	returned := int64(len(out) + len(tail))
+	got := r.GetStatistics()["samplesOut"]
+	if got != returned {
+		t.Errorf("samplesOut statistic = %d, want %d (Process %d + Flush %d): cubic Flush skips samplesOut accounting",
+			got, returned, len(out), len(tail))
+	}
+}
+
+// Process must return a caller-owned slice even when it emits nothing. During
+// the priming window Process yields zero output, and the empty slice it returns
+// must not alias the internal output buffer: a caller that appends to a
+// zero-length-but-nonzero-capacity result would otherwise have its data
+// silently overwritten by the next Process call reusing that buffer. This pins
+// the owned-empty-slice contract that the sibling stages already satisfy by
+// returning fresh []F{} literals (issue #51).
+func TestCubicStage_ProcessEmptyResultIsOwned(t *testing.T) {
+	c := NewCubicStage[float64](48000.0 / 44100.0)
+
+	// One sample leaves the stage partially primed (primed=1 <
+	// cubicLatencySamples=2), so Process emits nothing.
+	empty, err := c.Process([]float64{1.0})
+	if err != nil {
+		t.Fatal(err)
+	}
+	if len(empty) != 0 {
+		t.Fatalf("expected empty output during priming, got %d samples", len(empty))
+	}
+
+	// A caller reasonably appends its own data to the returned slice.
+	const sentinel = 12345.0
+	owned := append(empty, sentinel) //nolint:gocritic // deliberately appending to the returned slice to prove it is owned
+
+	// A subsequent Process that emits output must not corrupt the caller's
+	// slice. The input is small enough that Process reuses the same internal
+	// buffer without reallocating, so an aliased empty result would be
+	// overwritten here.
+	if _, err := c.Process([]float64{2.0, 3.0, 4.0, 5.0}); err != nil {
+		t.Fatal(err)
+	}
+
+	if owned[0] != sentinel {
+		t.Errorf("caller's appended value was corrupted: got %v, want %v (Process returned an aliased empty slice)", owned[0], sentinel)
+	}
+}
+
+// Cubic chunked-vs-one-shot equivalence at engine level: the root pin test
+// (streaming_equivalence_test.go) covers QualityLow/Medium/High only, not
+// QualityQuick's cubic path. Feeding the same signal in small chunks versus
+// one shot, with a single Flush at the end, must be bit-exact.
+func TestCubicStage_ChunkedEquivalence(t *testing.T) {
+	const n = 44100
+	in := make([]float64, n)
+	for i := range in {
+		in[i] = 0.5 * math.Sin(2*math.Pi*997*float64(i)/44100)
+	}
+
+	oneShot, err := NewResampler[float64](44100, 48000, QualityQuick)
+	if err != nil {
+		t.Fatal(err)
+	}
+	ref, err := oneShot.Process(slices.Clone(in))
+	if err != nil {
+		t.Fatal(err)
+	}
+	refTail, err := oneShot.Flush()
+	if err != nil {
+		t.Fatal(err)
+	}
+	ref = append(ref, refTail...)
+
+	plans := [][]int{
+		{1, 7, 13, 470, 4096},
+		{31, 331, 997},
+	}
+	for pi, plan := range plans {
+		chunked, err := NewResampler[float64](44100, 48000, QualityQuick)
+		if err != nil {
+			t.Fatal(err)
+		}
+		var got []float64
+		pos := 0
+		planIdx := 0
+		for pos < n {
+			size := plan[planIdx%len(plan)]
+			planIdx++
+			if pos+size > n {
+				size = n - pos
+			}
+			out, err := chunked.Process(slices.Clone(in[pos : pos+size]))
+			if err != nil {
+				t.Fatalf("plan %d: Process: %v", pi, err)
+			}
+			got = append(got, out...)
+			pos += size
+		}
+		tail, err := chunked.Flush()
+		if err != nil {
+			t.Fatalf("plan %d: Flush: %v", pi, err)
+		}
+		got = append(got, tail...)
+
+		if len(got) != len(ref) {
+			t.Fatalf("plan %d: length %d != one-shot %d", pi, len(got), len(ref))
+		}
+		for i := range got {
+			if got[i] != ref[i] {
+				t.Fatalf("plan %d: sample %d differs: %g != %g", pi, i, got[i], ref[i])
+			}
+		}
+	}
+}
diff --git a/internal/engine/debug_latency_test.go b/internal/engine/debug_latency_test.go
index 1c2b308..20c24f0 100644
--- a/internal/engine/debug_latency_test.go
+++ b/internal/engine/debug_latency_test.go
@@ -4,79 +4,57 @@
 package engine
 
 import (
-	"fmt"
+	"math"
 	"testing"
 )
 
-func TestDebugLatency(t *testing.T) {
-	resampler, err := NewResampler[float64](44100, 48000, QualityHigh)
+// TestLatencyReporting verifies Latency() (the startup deficit reported so
+// callers can prime a fixed-size output FIFO) against a deficit measured
+// independently of that accessor.
+//
+// The DFT pre-stage does valid-mode convolution (see processZeroCopy): it
+// never zero-pads history, so it withholds output entirely (returns an
+// empty slice) until enough real input has accumulated, rather than
+// emitting spurious zero-valued samples. That means a single one-shot
+// Process() call over a large chunk returns fewer output samples than an
+// ideal, delay-free resampler would (n*ratio): the shortfall is exactly the
+// output the filter is withholding while its delay lines prime, which
+// Flush() would later release. That shortfall is this test's independent
+// measurement of the startup deficit.
+func TestLatencyReporting(t *testing.T) {
+	r, err := NewResampler[float64](44100, 48000, QualityHigh)
 	if err != nil {
 		t.Fatal(err)
 	}
 
-	// Check stage parameters
-	fmt.Printf("\n=== Resampler Configuration ===\n")
-	fmt.Printf("Ratio: %.6f\n", resampler.ratio)
-
-	if resampler.preStage != nil {
-		fmt.Printf("\nDFT Pre-Stage:\n")
-		fmt.Printf("  Factor: %d\n", resampler.preStage.factor)
-		fmt.Printf("  Taps per phase: %d\n", resampler.preStage.tapsPerPhase)
-		fmt.Printf("  Filter latency: %d samples (input domain)\n", resampler.preStage.tapsPerPhase/2)
-	}
-
-	if resampler.polyphaseStage != nil {
-		fmt.Printf("\nPolyphase Stage:\n")
-		fmt.Printf("  Num phases: %d\n", resampler.polyphaseStage.numPhases)
-		fmt.Printf("  Taps per phase: %d\n", resampler.polyphaseStage.tapsPerPhase)
-		fmt.Printf("  Step: %d\n", resampler.polyphaseStage.step)
-		fmt.Printf("  Filter latency: %d samples (intermediate domain)\n", resampler.polyphaseStage.tapsPerPhase/2)
-	}
-
-	// Calculate total latency
-	// DFT stage latency (in input samples) + polyphase latency (in intermediate samples, converted to output)
-	dftLatency := 0
-	if resampler.preStage != nil && resampler.preStage.factor > 1 {
-		dftLatency = resampler.preStage.tapsPerPhase / 2
+	t.Logf("Resampler ratio: %.6f", r.ratio)
+	if r.preStage != nil {
+		t.Logf("DFT pre-stage: factor=%d tapsPerPhase=%d", r.preStage.factor, r.preStage.tapsPerPhase)
 	}
-
-	polyLatency := 0
-	if resampler.polyphaseStage != nil {
-		// Polyphase latency in intermediate samples
-		polyLatencyIntermediate := resampler.polyphaseStage.tapsPerPhase / 2
-		// Convert to output samples (approximately)
-		polyLatency = int(int64(polyLatencyIntermediate) * int64(resampler.polyphaseStage.numPhases) / resampler.polyphaseStage.step)
+	if r.polyphaseStage != nil {
+		t.Logf("Polyphase stage: numPhases=%d tapsPerPhase=%d step=%d",
+			r.polyphaseStage.numPhases, r.polyphaseStage.tapsPerPhase, r.polyphaseStage.step)
 	}
 
-	fmt.Printf("\nEstimated total latency: %d + %d ≈ %d output samples\n",
-		dftLatency*2, polyLatency, dftLatency*2+polyLatency)
-
-	// Test with small input to see when output appears
-	fmt.Printf("\n=== Testing Output Timing ===\n")
-	input := make([]float64, 100)
+	const n = 8192
+	input := make([]float64, n)
 	for i := range input {
-		input[i] = 1.0 // DC signal
+		input[i] = math.Cos(2 * math.Pi * 997 * float64(i) / 44100)
 	}
 
-	output, _ := resampler.Process(input)
-
-	// Count leading zeros
-	leadingZeros := 0
-	for _, v := range output {
-		if v == 0 {
-			leadingZeros++
-		} else {
-			break
-		}
+	output, err := r.Process(input)
+	if err != nil {
+		t.Fatalf("Process failed: %v", err)
 	}
 
-	fmt.Printf("Input: 100 samples\n")
-	fmt.Printf("Output: %d samples\n", len(output))
-	fmt.Printf("Leading zeros: %d\n", leadingZeros)
+	ideal := float64(n) * r.ratio
+	measuredDeficit := int(math.Round(ideal)) - len(output)
+
+	latency := r.Latency()
+	t.Logf("Input: %d samples, ideal output: %.1f, actual Process() output: %d, measured deficit: %d, Latency(): %d",
+		n, ideal, len(output), measuredDeficit, latency)
 
-	// Show first non-zero values
-	fmt.Printf("\nFirst 10 output values:\n")
-	for i := 0; i < 10 && i < len(output); i++ {
-		fmt.Printf("  [%d] %.10f\n", i, output[i])
+	if diff := latency - measuredDeficit; diff < -2 || diff > 2 {
+		t.Errorf("Latency() = %d, measured startup deficit = %d, want within 2 samples", latency, measuredDeficit)
 	}
 }
diff --git a/internal/engine/dft_stage.go b/internal/engine/dft_stage.go
index 3304f08..6fe2a8b 100644
--- a/internal/engine/dft_stage.go
+++ b/internal/engine/dft_stage.go
@@ -213,11 +213,10 @@ func (s *DFTStage[F]) Process(input []F) ([]F, error) {
 	if err != nil || len(output) == 0 {
 		return output, err
 	}
-	if s.factor == 1 {
-		return output, nil
-	}
-	// Return a copy to prevent caller's slice from being corrupted
-	// if they call Process() or Flush() again (which reuses s.outputBuf)
+	// processZeroCopy may return a slice that aliases the input (at factor==1,
+	// a passthrough) or the internal s.outputBuf; Process guarantees owned
+	// memory that stays valid across subsequent Process/Flush calls, so copy
+	// unconditionally.
 	result := make([]F, len(output))
 	copy(result, output)
 	return result, nil
@@ -227,13 +226,11 @@ func (s *DFTStage[F]) Process(input []F) ([]F, error) {
 // This is the simple path for small inputs.
 // Optimized with half-band detection, ConvolveValidMulti and Interleave2.
 func (s *DFTStage[F]) processChunk(history, output []F, numInputProcessable, factor, tapsPerPhase int) {
-	// Ensure phase buffers are large enough
+	// Ensure phase buffers are large enough. growStableLen adds growth slack
+	// when it must reallocate, matching the buffer policy used elsewhere in
+	// this file so one-sample chunk jitter does not trigger repeated reallocs.
 	for phase := range factor {
-		if cap(s.phaseBufs[phase]) < numInputProcessable {
-			s.phaseBufs[phase] = make([]F, numInputProcessable)
-		} else {
-			s.phaseBufs[phase] = s.phaseBufs[phase][:numInputProcessable]
-		}
+		s.phaseBufs[phase] = growStableLen(s.phaseBufs[phase], numInputProcessable)
 	}
 
 	historySlice := history[:numInputProcessable+tapsPerPhase-1]
@@ -343,9 +340,20 @@ func (s *DFTStage[F]) Flush() ([]F, error) {
 		return []F{}, nil
 	}
 
-	// Pad with zeros to flush pipeline
-	zeros := make([]F, s.tapsPerPhase)
-	return s.Process(zeros)
+	// Process retains exactly tapsPerPhase-1 history samples, so that many
+	// padding zeros advance the delay line past the last real sample without
+	// producing an extra all-zero output window.
+	zeros := make([]F, s.tapsPerPhase-1)
+	out, err := s.Process(zeros)
+	// Flush is terminal: after draining, return the stage to its fresh state.
+	// Otherwise the tapsPerPhase-1 padding zeros stay in the delay line, so the
+	// len(history)==0 guard never fires again: a second Flush keeps emitting
+	// all-zero windows and a post-flush Process convolves new audio against
+	// leftover zeros instead of starting a clean stream. Reset() is
+	// the authoritative fresh-state definition; calling it keeps Flush aligned
+	// with it automatically.
+	s.Reset()
+	return out, err
 }
 
 // Reset clears internal state.
@@ -560,13 +568,11 @@ func (s *DFTDecimationStage[F]) Process(input []F) ([]F, error) {
 	if err != nil || len(output) == 0 {
 		return output, err
 	}
-	if s.factor == 1 {
-		return output, nil
-	}
-	// IMPORTANT: Return a COPY of the output, not a slice of the internal buffer.
-	// Returning s.outputBuf directly would cause buffer corruption on the next
-	// Process() call, as the caller's slice would share the same backing array.
-	// This was the cause of TestResampler_BufferIntegrity failures for 96→48.
+	// processZeroCopy may return a slice that aliases the input (at factor==1,
+	// a passthrough) or the internal s.outputBuf; Process guarantees owned
+	// memory that stays valid across subsequent Process/Flush calls, so copy
+	// unconditionally. Returning s.outputBuf directly caused the
+	// TestResampler_BufferIntegrity failures for 96 to 48.
 	result := make([]F, len(output))
 	copy(result, output)
 	return result, nil
@@ -578,9 +584,20 @@ func (s *DFTDecimationStage[F]) Flush() ([]F, error) {
 		return []F{}, nil
 	}
 
-	// Pad with zeros to flush pipeline
-	zeros := make([]F, s.numTaps)
-	return s.Process(zeros)
+	// Process retains exactly numTaps-1 history samples, so that many padding
+	// zeros advance the delay line past the last real sample without producing
+	// an extra all-zero output window.
+	zeros := make([]F, s.numTaps-1)
+	out, err := s.Process(zeros)
+	// Flush is terminal: after draining, return the stage to its fresh state.
+	// Otherwise the numTaps-1 padding zeros stay in the delay line, so the
+	// len(history)==0 guard never fires again: a second Flush keeps emitting
+	// all-zero windows and a post-flush Process convolves new audio against
+	// leftover zeros instead of starting a clean stream. Reset() is
+	// the authoritative fresh-state definition (history and decimPhase); calling
+	// it keeps Flush aligned with it automatically.
+	s.Reset()
+	return out, err
 }
 
 // Reset clears internal state.
diff --git a/internal/engine/phase_wrap_measure_test.go b/internal/engine/phase_wrap_measure_test.go
new file mode 100644
index 0000000..06a46c8
--- /dev/null
+++ b/internal/engine/phase_wrap_measure_test.go
@@ -0,0 +1,370 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package engine
+
+import (
+	"math"
+	"math/cmplx"
+	"testing"
+)
+
+// Phase-wrap coefficient indexing investigation (measurement first).
+//
+// NewPolyphaseStage builds the cubic sub-phase interpolation banks (B/C/D) by
+// sampling the prototype at the current phase and its neighbours (getCoeff in
+// polyphase_stage.go). Two reviews argued that the original code, which wrapped
+// the neighbour phase WITHIN the same tap (phase % numPhases), picked a
+// prototype sample numPhases-1 positions away at each phase boundary instead of
+// the adjacent one. The mathematically adjacent sample is the neighbour in the
+// FLAT prototype: the right neighbour of (tap t, phase L-1) is prototype[t*L +
+// L] (phase 0 of tap t+1), and the left neighbour of (tap t, phase 0) is
+// prototype[t*L - 1] (phase L-1 of tap t-1); out-of-range positions clamp to
+// 0.0.
+//
+// This test measured both indexings on the canonical 997 Hz full-scale sine and
+// drove the decision to switch getCoeff to flat indexing. It is retained as the
+// record of that measurement and as a regression guard: if getCoeff is ever
+// reverted to the wrap, the bank-equality check and the active-ratio THD+N
+// assertion below both fail.
+//
+// KEY FINDING: the fixed-point sub-phase accumulator only produces a non-zero
+// fractional part x (which is what the B/C/D banks interpolate over) when the
+// resampling ratio is NOT an exact rational the phase count L divides cleanly.
+// For 44100<->48000 the polyphase stage ratio is exactly 80/147, so the step is
+// an integer number of phase units, x is identically 0 for every output sample,
+// and the B/C/D banks (hence the disputed wrap) are never consulted. That is why
+// the wrap error was invisible to the entire existing regression suite, whose
+// ratios (44100<->48000, 48000->32000, 48000->96000) are all exact-rational or
+// integer/DFT-only. The measurement therefore uses TWO probes: the degenerate
+// exact ratio (to document the x==0 masking) and an active-interpolation ratio
+// (to actually reach the code under test).
+//
+// The 44100 -> N QualityHigh pipeline is a 2x DFT pre-stage (44100 -> 88200)
+// followed by the polyphase stage, so the measurement drives that stage in
+// isolation with a clean full-scale sine at the intermediate 88200 Hz rate: a
+// pure tone is the ideal probe for coefficient-interpolation distortion, and
+// isolating the stage keeps everything identical between runs except the
+// boundary indexing.
+
+const (
+	phaseWrapPreUpsample = 2
+	phaseWrapTestFreq    = 997.0
+	phaseWrapAmplitude   = 1.0    // full-scale (0 dBFS)
+	phaseWrapNumSamples  = 131072 // input samples at the intermediate 88200 Hz rate
+	phaseWrapFFTSize     = 32768
+	phaseWrapFFTOffset   = 8192 // skip the filter start-up transient
+	phaseWrapFundGuard   = 12   // bins around the fundamental treated as signal
+	phaseWrapDCGuard     = 3    // low bins excluded as window DC leakage
+
+	// Regression bounds on the active-interpolation probe (44100 -> 64000).
+	// MEASURED on this machine (float64, values below): the production flat path
+	// yields about -140.7 dB THD+N and the old wrap about -54.5 dB, an improvement
+	// of about 86 dB. The bounds leave generous headroom so the guard fires only on
+	// a genuine regression (e.g. reverting getCoeff to the wrap), not on minor
+	// numeric drift.
+	phaseWrapFlatCeilingDB    = -120.0 // production (flat) THD+N must be at least this clean
+	phaseWrapImprovementFloor = 40.0   // flat must beat the old wrap by at least this much
+	phaseWrapPlausibleFloorDB = -160.0 // below this, the harness is broken, not the library
+	phaseWrapPlausibleCeilDB  = -40.0  // above this for flat, the harness is broken
+)
+
+// phaseWrapResult holds the outcome of a single ratio measurement.
+type phaseWrapResult struct {
+	flatTHDN   float64 // THD+N of the production (flat) path
+	wrapTHDN   float64 // THD+N of the old wrap path
+	sampleDiff int     // number of output samples that differ between the two
+	total      int     // total compared samples
+	numPhases  int
+}
+
+// buildPhaseInterpBanks reproduces the cubic interpolation bank construction
+// from NewPolyphaseStage exactly, parameterised only by the boundary indexing.
+// With flat=true it reproduces the production (fixed) indexing; with flat=false
+// it reproduces the original wrap. Everything else, the Catmull-Rom cubic math
+// and the reversed tap storage, is identical, so any output difference is
+// attributable solely to the boundary indexing.
+func buildPhaseInterpBanks(fb *polyphaseFilter, numPhases int, flat bool) (a, b, c, d [][]float64) {
+	tapsPerPhase := fb.tapsPerPhase
+	coeffs := fb.coeffs
+
+	getCoeff := func(phase, tap int) float64 {
+		var idx int
+		if flat {
+			// Flat-prototype-spill: neighbour of a phase boundary is the adjacent
+			// prototype sample, crossing into the next/previous tap. Matches the
+			// production getCoeff after the fix.
+			idx = tap*numPhases + phase
+		} else {
+			// Original wrap: keep the neighbour within the same tap.
+			wrappedPhase := phase % numPhases
+			if wrappedPhase < 0 {
+				wrappedPhase += numPhases
+			}
+			idx = tap*numPhases + wrappedPhase
+		}
+		if idx < 0 || idx >= len(coeffs) {
+			return 0.0
+		}
+		return coeffs[idx]
+	}
+
+	a = make([][]float64, numPhases)
+	b = make([][]float64, numPhases)
+	c = make([][]float64, numPhases)
+	d = make([][]float64, numPhases)
+
+	for phase := range numPhases {
+		a[phase] = make([]float64, tapsPerPhase)
+		b[phase] = make([]float64, tapsPerPhase)
+		c[phase] = make([]float64, tapsPerPhase)
+		d[phase] = make([]float64, tapsPerPhase)
+
+		for tap := range tapsPerPhase {
+			f0 := getCoeff(phase, tap)
+			f1 := getCoeff(phase+1, tap)
+			fm1 := getCoeff(phase-1, tap)
+			f2 := getCoeff(phase+cubicPhaseOffset, tap)
+
+			av := f0
+			cv := cubicCenterCoeff*(f1+fm1) - f0
+			dv := (1.0 / cubicDivisor) * (f2 - f1 + fm1 - f0 - cubicCMultiplier*cv)
+			bv := f1 - f0 - dv - cv
+
+			revTap := tapsPerPhase - 1 - tap
+			a[phase][revTap] = av
+			b[phase][revTap] = bv
+			c[phase][revTap] = cv
+			d[phase][revTap] = dv
+		}
+	}
+	return a, b, c, d
+}
+
+// processAllPhaseWrap runs Process followed by the terminal Flush and returns
+// the full output stream.
+func processAllPhaseWrap(t *testing.T, s *PolyphaseStage[float64], input []float64) []float64 {
+	t.Helper()
+	out, err := s.Process(input)
+	if err != nil {
+		t.Fatalf("Process failed: %v", err)
+	}
+	flush, err := s.Flush()
+	if err != nil {
+		t.Fatalf("Flush failed: %v", err)
+	}
+	return append(out, flush...)
+}
+
+// blackmanHarris7 returns a 7-term Blackman-Harris window of length n. Its peak
+// sidelobe is about -180 dB, far below the THD+N floor of interest, so spectral
+// leakage from the non-bin-aligned 997 Hz tone does not mask the measurement.
+func blackmanHarris7(n int) []float64 {
+	const (
+		a0 = 0.27105140069342
+		a1 = 0.43329793923448
+		a2 = 0.21812299954311
+		a3 = 0.06592544638803
+		a4 = 0.01081174209837
+		a5 = 0.00077658482522
+		a6 = 0.00001388721735
+	)
+	w := make([]float64, n)
+	denom := float64(n - 1)
+	for i := range w {
+		t := 2.0 * math.Pi * float64(i) / denom
+		w[i] = a0 - a1*math.Cos(t) + a2*math.Cos(2*t) - a3*math.Cos(3*t) +
+			a4*math.Cos(4*t) - a5*math.Cos(5*t) + a6*math.Cos(6*t)
+	}
+	return w
+}
+
+// measureTHDNPhaseWrap computes THD+N in dB for a resampled tone. It notches the
+// fundamental (and its window main lobe) plus the DC region, then ratios all
+// remaining spectral power (harmonics, images, noise) against the fundamental
+// power. A more negative figure is better.
+func measureTHDNPhaseWrap(t *testing.T, output []float64, testFreq, sampleRate float64) float64 {
+	t.Helper()
+	if len(output) < phaseWrapFFTOffset+phaseWrapFFTSize {
+		t.Fatalf("output too short for THD+N: have %d, need %d", len(output), phaseWrapFFTOffset+phaseWrapFFTSize)
+	}
+
+	win := blackmanHarris7(phaseWrapFFTSize)
+	fftIn := make([]complex128, phaseWrapFFTSize)
+	for i := range phaseWrapFFTSize {
+		fftIn[i] = complex(output[phaseWrapFFTOffset+i]*win[i], 0)
+	}
+	fftOut := fft(fftIn)
+
+	half := phaseWrapFFTSize / 2
+	fundBin := int(math.Round(testFreq / sampleRate * float64(phaseWrapFFTSize)))
+
+	var signalPower, totalPower float64
+	for bin := phaseWrapDCGuard + 1; bin < half; bin++ {
+		mag := cmplx.Abs(fftOut[bin])
+		p := mag * mag
+		totalPower += p
+		if bin >= fundBin-phaseWrapFundGuard && bin <= fundBin+phaseWrapFundGuard {
+			signalPower += p
+		}
+	}
+
+	noiseDistPower := totalPower - signalPower
+	if signalPower <= 0 {
+		t.Fatalf("no signal power measured (fundBin=%d)", fundBin)
+	}
+	return 10.0 * math.Log10(noiseDistPower/signalPower+1e-30)
+}
+
+// measurePhaseWrap resamples a full-scale 997 Hz sine through the polyphase
+// stage for inputRate -> outputRate (QualityHigh) with both the production flat
+// indexing and the old wrap, and returns the THD+N of each plus how many output
+// samples differ. It also asserts that the production stage's banks match the
+// flat reconstruction, proving the fix is in place.
+func measurePhaseWrap(t *testing.T, inputRate, outputRate float64) phaseWrapResult {
+	t.Helper()
+
+	intermediate := inputRate * phaseWrapPreUpsample // rate the polyphase stage sees
+	polyphaseRatio := outputRate / intermediate
+	totalIORatio := inputRate / outputRate
+	const hasPreStage = true // 44100 -> N here is non-integer upsampling: 2x DFT pre-stage
+
+	numPhases, _ := findRationalApprox(polyphaseRatio)
+
+	fb, err := designPolyphaseFilter(numPhases, polyphaseRatio, totalIORatio, hasPreStage, QualityHigh)
+	if err != nil {
+		t.Fatalf("designPolyphaseFilter failed: %v", err)
+	}
+
+	flatA, flatB, flatC, flatD := buildPhaseInterpBanks(fb, numPhases, true)
+	wrapA, wrapB, wrapC, wrapD := buildPhaseInterpBanks(fb, numPhases, false)
+
+	// Production stage uses the fixed flat indexing; verify its banks match the
+	// flat reconstruction. If getCoeff is ever reverted to the wrap this fails.
+	flatStage, err := NewPolyphaseStage[float64](polyphaseRatio, totalIORatio, hasPreStage, QualityHigh)
+	if err != nil {
+		t.Fatalf("NewPolyphaseStage (flat) failed: %v", err)
+	}
+	assertBanksEqual(t, "polyCoeffs", flatStage.polyCoeffs, flatA)
+	assertBanksEqual(t, "polyCoeffsB", flatStage.polyCoeffsB, flatB)
+	assertBanksEqual(t, "polyCoeffsC", flatStage.polyCoeffsC, flatC)
+	assertBanksEqual(t, "polyCoeffsD", flatStage.polyCoeffsD, flatD)
+
+	// Second stage identical to the first except its banks are overwritten with
+	// the old wrap variant.
+	wrapStage, err := NewPolyphaseStage[float64](polyphaseRatio, totalIORatio, hasPreStage, QualityHigh)
+	if err != nil {
+		t.Fatalf("NewPolyphaseStage (wrap) failed: %v", err)
+	}
+	wrapStage.polyCoeffs = wrapA
+	wrapStage.polyCoeffsB = wrapB
+	wrapStage.polyCoeffsC = wrapC
+	wrapStage.polyCoeffsD = wrapD
+
+	// Full-scale 997 Hz sine sampled at the intermediate rate.
+	input := make([]float64, phaseWrapNumSamples)
+	for i := range input {
+		input[i] = phaseWrapAmplitude * math.Sin(2.0*math.Pi*phaseWrapTestFreq*float64(i)/intermediate)
+	}
+
+	flatOut := processAllPhaseWrap(t, flatStage, input)
+	wrapOut := processAllPhaseWrap(t, wrapStage, input)
+
+	n := min(len(flatOut), len(wrapOut))
+	sampleDiff := 0
+	for i := range n {
+		if flatOut[i] != wrapOut[i] {
+			sampleDiff++
+		}
+	}
+
+	return phaseWrapResult{
+		flatTHDN:   measureTHDNPhaseWrap(t, flatOut, phaseWrapTestFreq, outputRate),
+		wrapTHDN:   measureTHDNPhaseWrap(t, wrapOut, phaseWrapTestFreq, outputRate),
+		sampleDiff: sampleDiff,
+		total:      n,
+		numPhases:  numPhases,
+	}
+}
+
+// TestPhaseWrapCoeffIndexing_Measure records the measurement that drove the
+// getCoeff wrap -> flat fix and guards against a regression.
+//
+// MEASURED RESULT (this machine, float64):
+//
+//	Probe 1, 44100->48000 (exact 80/147, sub-phase x == 0):
+//	  wrap THD+N = -140.72 dB, flat THD+N = -140.72 dB, improvement = 0.00 dB
+//	  0 / 71332 output samples differ (bit-identical: the B/C/D banks are never
+//	  consulted, so the wrap error cannot manifest at this ratio).
+//
+//	Probe 2, 44100->64000 (active sub-phase interpolation, L = 201):
+//	  wrap THD+N = -54.46 dB, flat THD+N = -140.72 dB, improvement = 86.26 dB
+//	  1486 / 95109 output samples differ.
+//
+// The flat path (~-140.7 dB) matches soxr's HQ-class THD+N; the wrap path
+// (-54 dB) is wildly worse, confirming the wrap was a genuine defect rather than
+// a harness artifact. The decision to apply the fix follows: flat improves
+// THD+N by far more than the 1 dB plan threshold wherever the disputed code is
+// reachable. The 44100->48000 probe reads 0 dB only because that ratio's
+// sub-phase interpolation is inert (x == 0), which is exactly why the bug hid.
+func TestPhaseWrapCoeffIndexing_Measure(t *testing.T) {
+	// Probe 1: degenerate exact ratio. Documents the x == 0 masking.
+	deg := measurePhaseWrap(t, 44100, 48000)
+	if deg.numPhases != 80 {
+		t.Fatalf("expected 80 phases for 48000/88200, got %d", deg.numPhases)
+	}
+	t.Logf("44100->48000 (exact 80/147, sub-phase x==0): wrap THD+N=%.2f dB, flat THD+N=%.2f dB, improvement=%.2f dB; %d/%d samples differ",
+		deg.wrapTHDN, deg.flatTHDN, deg.wrapTHDN-deg.flatTHDN, deg.sampleDiff, deg.total)
+	if deg.sampleDiff != 0 {
+		t.Errorf("expected bit-identical wrap/flat output at exact ratio 44100->48000 (x==0), got %d differing samples", deg.sampleDiff)
+	}
+
+	// Probe 2: active sub-phase interpolation. Reaches the disputed code.
+	act := measurePhaseWrap(t, 44100, 64000)
+	improvement := act.wrapTHDN - act.flatTHDN
+	t.Logf("44100->64000 (active sub-phase interp, L=%d): wrap THD+N=%.2f dB, flat THD+N=%.2f dB, improvement=%.2f dB; %d/%d samples differ",
+		act.numPhases, act.wrapTHDN, act.flatTHDN, improvement, act.sampleDiff, act.total)
+
+	if act.sampleDiff == 0 {
+		t.Fatalf("expected wrap/flat output to differ at active ratio 44100->64000; harness not exercising the code under test")
+	}
+
+	// Sanity: the production flat THD+N must be in a plausible HQ range. A
+	// wildly-off figure means the harness is broken, not that the library is that
+	// good or bad.
+	if act.flatTHDN < phaseWrapPlausibleFloorDB || act.flatTHDN > phaseWrapPlausibleCeilDB {
+		t.Fatalf("flat THD+N %.2f dB implausible; measurement harness likely broken", act.flatTHDN)
+	}
+
+	// Regression: the production (flat) path must stay clean.
+	if act.flatTHDN > phaseWrapFlatCeilingDB {
+		t.Errorf("flat THD+N regression: got %.2f dB, want <= %.2f dB", act.flatTHDN, phaseWrapFlatCeilingDB)
+	}
+
+	// Regression: the fix must keep beating the old wrap by a wide margin. If this
+	// shrinks, getCoeff may have been reverted to the wrap.
+	if improvement < phaseWrapImprovementFloor {
+		t.Errorf("flat improves THD+N by only %.2f dB over the old wrap, want >= %.2f dB; "+
+			"getCoeff may have regressed to the phase wrap", improvement, phaseWrapImprovementFloor)
+	}
+}
+
+// assertBanksEqual fails the test if two coefficient banks differ, proving the
+// production stage's banks match the flat reconstruction exactly.
+func assertBanksEqual(t *testing.T, name string, got, want [][]float64) {
+	t.Helper()
+	if len(got) != len(want) {
+		t.Fatalf("%s: phase count mismatch: got %d want %d", name, len(got), len(want))
+	}
+	for phase := range want {
+		if len(got[phase]) != len(want[phase]) {
+			t.Fatalf("%s: phase %d tap count mismatch: got %d want %d", name, phase, len(got[phase]), len(want[phase]))
+		}
+		for tap := range want[phase] {
+			if got[phase][tap] != want[phase][tap] {
+				t.Fatalf("%s: mismatch at phase %d tap %d: got %v want %v",
+					name, phase, tap, got[phase][tap], want[phase][tap])
+			}
+		}
+	}
+}
diff --git a/internal/engine/polyphase_stage.go b/internal/engine/polyphase_stage.go
index 834457b..497d0e9 100644
--- a/internal/engine/polyphase_stage.go
+++ b/internal/engine/polyphase_stage.go
@@ -101,14 +101,28 @@ func NewPolyphaseStage[F simdops.Float](ratio, totalIORatio float64, hasPreStage
 	phaseFracScale := float64(int64(1) << phaseFracBits)
 	step := int64(math.Round((1.0 / ratio) * float64(numPhases) * phaseFracScale))
 
-	// Helper function to get prototype coefficient with wrap-around for interpolation
+	// getCoeff samples the prototype for cubic sub-phase interpolation of the
+	// polyphase banks. filterBank.coeffs is the flat prototype
+	// (coeffs[tap*numPhases + phase] == prototype[tap*numPhases + phase]), so the
+	// coefficient adjacent to a phase boundary is the neighbour in this FLAT
+	// array: the right neighbour of (tap t, phase L-1) is coeffs[t*L + L], i.e.
+	// phase 0 of tap t+1, and the left neighbour of (tap t, phase 0) is
+	// coeffs[t*L - 1], i.e. phase L-1 of tap t-1. Boundary policy: positions
+	// before the first or after the last prototype sample have no data, so they
+	// clamp to 0.0 (the natural tails of a finite impulse response).
+	//
+	// This must NOT wrap phase within the same tap (phase % numPhases). Wrapping
+	// picks a prototype sample numPhases-1 positions away at each boundary, which
+	// injects a large discontinuity into the interpolated coefficient. That error
+	// is invisible for exact-rational ratios (e.g. 44100<->48000 == 80/147, where
+	// the fixed-point sub-phase x is identically 0 so the B/C/D banks are never
+	// consulted), but for ratios with active sub-phase interpolation it collapses
+	// THD+N by a large margin: the committed measurement (phase_wrap_measure_test.go)
+	// reads -54.46 dB wrapped versus -140.72 dB flat at 44100 -> 64000, an 86.26 dB
+	// improvement, with similar magnitude at other active-interpolation ratios
+	// measured during the investigation but not committed as tests.
 	getCoeff := func(phase, tap int) float64 {
-		// Wrap phase around for interpolation at boundaries
-		wrappedPhase := phase % numPhases
-		if wrappedPhase < 0 {
-			wrappedPhase += numPhases
-		}
-		idx := tap*numPhases + wrappedPhase
+		idx := tap*numPhases + phase
 		if idx < 0 || idx >= len(filterBank.coeffs) {
 			return 0.0
 		}
@@ -295,16 +309,19 @@ func (s *PolyphaseStage[F]) processZeroCopy(input []F) ([]F, error) { //nolint:u
 	// Trim output to actual size produced
 	output := s.outputBuf[:outIdx]
 
-	// Consume processed samples from history
-	consumed := int(at>>phaseFracBits) / numPhases
-	if consumed > 0 && consumed <= histLen {
+	// Consume processed samples from history. The accumulator can overshoot
+	// limit by up to one step at severe downsampling ratios, pointing past the
+	// fully available input positions; cap at numIn so the delay line always
+	// retains tapsPerPhase-1 samples and the rebase below matches the trim.
+	consumed := min(int((at>>phaseFracBits)/numPhases64), numIn)
+	if consumed > 0 {
 		copy(s.history, s.history[consumed:])
 		s.history = s.history[:histLen-consumed]
 	}
 
-	// Save remainder for next call
-	// Keep the fractional part within one input sample
-	s.at = at - int64(consumed*numPhases)<<phaseFracBits
+	// Save remainder for next call; may legitimately point past the retained
+	// history when the ratio is severe, the next call's input covers it.
+	s.at = at - int64(consumed)*numPhases64<<phaseFracBits
 
 	s.samplesOut += int64(len(output))
 
@@ -333,14 +350,21 @@ func (s *PolyphaseStage[F]) Flush() ([]F, error) {
 		return []F{}, nil
 	}
 
-	// Pad with tapsPerPhase zeros to drain the polyphase delay line, matching
-	// the sibling DFTStage.Flush. historyBufferMultiplier is a buffer
-	// pre-allocation constant, not a flush-padding amount: padding
-	// tapsPerPhase*historyBufferMultiplier zeros pushes an extra tapsPerPhase
-	// zeros through the filter, producing additional all-zero output windows
-	// (trailing silence) and a longer-than-canonical output length (issue #30).
-	zeros := make([]F, s.tapsPerPhase)
-	return s.Process(zeros)
+	// Process retains exactly tapsPerPhase-1 history samples, so that many
+	// padding zeros advance the delay line past the last real sample without
+	// producing an extra all-zero output window. historyBufferMultiplier is a
+	// buffer pre-allocation constant, not a flush-padding amount.
+	zeros := make([]F, s.tapsPerPhase-1)
+	out, err := s.Process(zeros)
+	// Flush is terminal: after draining, return the stage to its fresh state.
+	// Otherwise the tapsPerPhase-1 padding zeros stay in the delay line, so the
+	// len(history)==0 guard never fires again: a second Flush keeps emitting
+	// all-zero windows and a post-flush Process convolves new audio against
+	// leftover zeros instead of starting a clean stream. Reset() is
+	// the authoritative fresh-state definition (phase accumulator, history, and
+	// sample counters); calling it keeps Flush aligned with it automatically.
+	s.Reset()
+	return out, err
 }
 
 // Reset clears internal state.
diff --git a/internal/engine/resampler.go b/internal/engine/resampler.go
index 4879ec8..75b9042 100644
--- a/internal/engine/resampler.go
+++ b/internal/engine/resampler.go
@@ -49,8 +49,8 @@ type Resampler[F simdops.Float] struct {
 //   - For integer ratios: Uses only DFT stage
 //   - For non-integer ratios: Uses DFT pre-stage (2×) + polyphase stage
 func NewResampler[F simdops.Float](inputRate, outputRate float64, quality Quality) (*Resampler[F], error) {
-	if inputRate <= 0 || outputRate <= 0 {
-		return nil, fmt.Errorf("sample rates must be positive: input=%f, output=%f", inputRate, outputRate)
+	if !(inputRate > 0) || !(outputRate > 0) {
+		return nil, fmt.Errorf("sample rates must be positive finite numbers: input=%f, output=%f", inputRate, outputRate)
 	}
 
 	ratio := outputRate / inputRate
@@ -59,9 +59,11 @@ func NewResampler[F simdops.Float](inputRate, outputRate float64, quality Qualit
 	// Following SOXR's pattern: ratios between 1/256 and 256 are practical for audio.
 	// Extreme ratios can cause: (1) integer overflow in output size calculation,
 	// (2) memory exhaustion from attempting to allocate huge output buffers.
+	// The checks are expressed positively so a NaN ratio (e.g. from Inf/Inf)
+	// fails validation instead of silently passing every comparison.
 	const minRatio = 1.0 / 256.0 // 256x downsampling
 	const maxRatio = 256.0       // 256x upsampling
-	if ratio < minRatio || ratio > maxRatio {
+	if !(ratio >= minRatio && ratio <= maxRatio) {
 		return nil, fmt.Errorf("resampling ratio %.6f out of valid range [%.6f, %.0f]", ratio, minRatio, maxRatio)
 	}
 	ops := simdops.For[F]()
@@ -237,7 +239,9 @@ func (r *Resampler[F]) ProcessZeroCopy(input []F) ([]F, error) { //nolint:dupl /
 	r.samplesIn += int64(len(input))
 
 	if r.cubicStage != nil {
-		output, err := r.cubicStage.Process(input)
+		// Zero-copy path: the returned slice may alias the cubic stage's
+		// internal output buffer, matching the ProcessZeroCopy contract.
+		output, err := r.cubicStage.processZeroCopy(input)
 		if err != nil {
 			return nil, fmt.Errorf("cubic stage processing failed: %w", err)
 		}
@@ -273,9 +277,18 @@ func (r *Resampler[F]) ProcessZeroCopy(input []F) ([]F, error) { //nolint:dupl /
 
 // Flush returns any remaining buffered samples.
 func (r *Resampler[F]) Flush() ([]F, error) {
-	// QualityQuick cubic stage doesn't buffer
+	// QualityQuick cubic stage holds a cubicLatencySamples-sample tail; drain it.
+	// The samplesOut accounting at the bottom of this method is on the FIR
+	// path only, so this branch must account for its own emitted tail:
+	// otherwise GetStatistics undercounts by the tail length now that cubic
+	// Flush emits a real tail instead of always being empty.
 	if r.cubicStage != nil {
-		return r.cubicStage.Flush()
+		output, err := r.cubicStage.Flush()
+		if err != nil {
+			return nil, err
+		}
+		r.samplesOut += int64(len(output))
+		return output, nil
 	}
 
 	var output []F
@@ -352,6 +365,31 @@ func (r *Resampler[F]) GetStatistics() map[string]int64 {
 	}
 }
 
+// Latency returns the startup deficit in output samples: how many samples
+// the first Process calls withhold while the filter delay lines prime.
+// Real-time consumers should prime their output FIFO with this many samples
+// of silence to keep fixed-size callbacks fed.
+func (r *Resampler[F]) Latency() int {
+	if r.cubicStage != nil {
+		return int(math.Ceil(cubicLatencySamples * r.ratio))
+	}
+	deficitIn := 0.0
+	if r.preStage != nil && r.preStage.factor > 1 {
+		deficitIn += float64(r.preStage.tapsPerPhase - 1)
+	}
+	if r.decimationStage != nil {
+		deficitIn += float64(r.decimationStage.numTaps - 1)
+	}
+	if r.polyphaseStage != nil {
+		intermediateFactor := 1.0
+		if r.preStage != nil && r.preStage.factor > 1 {
+			intermediateFactor = float64(r.preStage.factor)
+		}
+		deficitIn += float64(r.polyphaseStage.tapsPerPhase-1) / intermediateFactor
+	}
+	return int(math.Ceil(deficitIn * r.ratio))
+}
+
 // isIntegerRatio checks if the ratio is an integer (within tolerance).
 func isIntegerRatio(ratio float64) bool {
 	const tolerance = 1e-9
diff --git a/internal/engine/reset_state_test.go b/internal/engine/reset_state_test.go
index 9b35e12..894ed7d 100644
--- a/internal/engine/reset_state_test.go
+++ b/internal/engine/reset_state_test.go
@@ -103,6 +103,7 @@ func TestResampler_Reset(t *testing.T) {
 		{"44100_to_48000", 44100, 48000},
 		{"44100_to_96000", 44100, 96000},
 		{"44100_to_88200", 44100, 88200},
+		{"48000_to_16000", 48000, 16000}, // integer decimation ratio, exercises decimationStage.Reset()
 	}
 
 	for _, tc := range testCases {
@@ -212,7 +213,10 @@ func TestInterpolationStages_Reset(t *testing.T) {
 // Multiple Reset Tests - Verify Reset() can be called multiple times
 // =============================================================================
 
-// TestDFTStage_MultipleResets verifies Reset() can be called repeatedly.
+// TestDFTStage_MultipleResets verifies Reset() can be called repeatedly, and
+// that identical input reproduces round 0's output bit-exactly after every
+// Reset() (a partially-cleared delay line would drift round to round even
+// though each individual output stays NaN/Inf-free).
 func TestDFTStage_MultipleResets(t *testing.T) {
 	stage, err := NewDFTStage[float64](2, QualityHigh)
 	require.NoError(t, err)
@@ -222,6 +226,7 @@ func TestDFTStage_MultipleResets(t *testing.T) {
 		input[i] = 1.0
 	}
 
+	var round0Output []float64
 	for round := range 5 {
 		// Process
 		output, err := stage.Process(input)
@@ -232,14 +237,29 @@ func TestDFTStage_MultipleResets(t *testing.T) {
 			assert.False(t, math.IsNaN(v), "Round %d: output[%d] is NaN", round, i)
 		}
 
+		if round == 0 {
+			round0Output = make([]float64, len(output))
+			copy(round0Output, output)
+		} else {
+			require.Len(t, output, len(round0Output), "Round %d: output length differs from round 0", round)
+			for i := range output {
+				// require stops at the first mismatch instead of flooding
+				// the log across an entire round's worth of samples.
+				require.InDelta(t, round0Output[i], output[i], 1e-15,
+					"Round %d: output[%d] differs from round 0 after Reset()", round, i)
+			}
+		}
+
 		// Reset
 		stage.Reset()
 	}
 
-	t.Log("DFT stage: 5 Reset() cycles completed successfully")
+	t.Log("DFT stage: 5 Reset() cycles reproduce round 0 output bit-exactly")
 }
 
-// TestResampler_MultipleResets verifies Resampler Reset() can be called repeatedly.
+// TestResampler_MultipleResets verifies Resampler Reset() can be called
+// repeatedly, and that identical input reproduces round 0's output
+// bit-exactly after every Reset().
 func TestResampler_MultipleResets(t *testing.T) {
 	resampler, err := NewResampler[float64](44100, 48000, QualityHigh)
 	require.NoError(t, err)
@@ -249,6 +269,7 @@ func TestResampler_MultipleResets(t *testing.T) {
 		input[i] = math.Sin(2.0 * math.Pi * float64(i) / 100)
 	}
 
+	var round0Output []float64
 	for round := range 5 {
 		// Process
 		output, err := resampler.Process(input)
@@ -259,11 +280,105 @@ func TestResampler_MultipleResets(t *testing.T) {
 			assert.False(t, math.IsNaN(v), "Round %d: output[%d] is NaN", round, i)
 		}
 
+		if round == 0 {
+			round0Output = make([]float64, len(output))
+			copy(round0Output, output)
+		} else {
+			require.Len(t, output, len(round0Output), "Round %d: output length differs from round 0", round)
+			for i := range output {
+				require.InDelta(t, round0Output[i], output[i], 1e-15,
+					"Round %d: output[%d] differs from round 0 after Reset()", round, i)
+			}
+		}
+
 		// Reset
 		resampler.Reset()
 	}
 
-	t.Log("Resampler: 5 Reset() cycles completed successfully")
+	t.Log("Resampler: 5 Reset() cycles reproduce round 0 output bit-exactly")
+}
+
+// TestResampler_ResetAfterFlush verifies that Reset() following a terminal
+// Flush() leaves the resampler equivalent to a fresh instance: processing
+// the same input again after Reset must match a fresh resampler's Process()
+// and Flush() bit-exactly. This guards the case a plain mid-stream Reset
+// test cannot: Flush() itself mutates state (draining delay lines, marking
+// the stream terminal), so Reset() has more to undo here than after a bare
+// Process().
+func TestResampler_ResetAfterFlush(t *testing.T) {
+	resampler, err := NewResampler[float64](44100, 48000, QualityHigh)
+	require.NoError(t, err)
+
+	input := make([]float64, 2000)
+	for i := range input {
+		input[i] = math.Sin(2.0 * math.Pi * 1000 * float64(i) / 44100)
+	}
+
+	_, err = resampler.Process(input)
+	require.NoError(t, err)
+	_, err = resampler.Flush()
+	require.NoError(t, err)
+
+	resampler.Reset()
+
+	output, err := resampler.Process(input)
+	require.NoError(t, err)
+	flush, err := resampler.Flush()
+	require.NoError(t, err)
+
+	freshResampler, err := NewResampler[float64](44100, 48000, QualityHigh)
+	require.NoError(t, err)
+	freshOutput, err := freshResampler.Process(input)
+	require.NoError(t, err)
+	freshFlush, err := freshResampler.Flush()
+	require.NoError(t, err)
+
+	require.Len(t, output, len(freshOutput), "Process length after Reset-following-Flush differs from fresh")
+	for i := range output {
+		require.InDelta(t, freshOutput[i], output[i], 1e-15,
+			"Process[%d] after Reset-following-Flush differs from fresh", i)
+	}
+	require.Len(t, flush, len(freshFlush), "Flush length after Reset-following-Flush differs from fresh")
+	for i := range flush {
+		require.InDelta(t, freshFlush[i], flush[i], 1e-15,
+			"Flush[%d] after Reset-following-Flush differs from fresh", i)
+	}
+
+	t.Log("Resampler: Reset() after terminal Flush() reproduces fresh-instance output bit-exactly")
+}
+
+// TestResampler_Reset_Float32 mirrors TestResampler_Reset for the float32
+// instantiation: Reset() must clear state correctly on the SIMD float32
+// code path too, not just float64.
+func TestResampler_Reset_Float32(t *testing.T) {
+	resampler, err := NewResampler[float32](44100, 48000, QualityHigh)
+	require.NoError(t, err)
+
+	input := make([]float32, 4000)
+	for i := range input {
+		input[i] = float32(math.Sin(2.0 * math.Pi * 1000 * float64(i) / 44100))
+	}
+	output1, err := resampler.Process(input)
+	require.NoError(t, err)
+	require.NotEmpty(t, output1)
+
+	resampler.Reset()
+
+	output2, err := resampler.Process(input)
+	require.NoError(t, err)
+
+	freshResampler, err := NewResampler[float32](44100, 48000, QualityHigh)
+	require.NoError(t, err)
+	outputFresh, err := freshResampler.Process(input)
+	require.NoError(t, err)
+
+	require.Len(t, output2, len(outputFresh), "Output length after reset should match fresh resampler")
+	for i := range output2 {
+		require.InDelta(t, outputFresh[i], output2[i], 1e-15,
+			"Output[%d] after reset differs from fresh resampler", i)
+	}
+
+	t.Log("Resampler[float32]: Reset() reproduces fresh-instance output bit-exactly")
 }
 
 // =============================================================================
diff --git a/internal/engine/severe_ratio_test.go b/internal/engine/severe_ratio_test.go
new file mode 100644
index 0000000..1fb5d74
--- /dev/null
+++ b/internal/engine/severe_ratio_test.go
@@ -0,0 +1,82 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package engine
+
+import (
+	"testing"
+)
+
+// Severe non-integer downsampling: the fixed-point accumulator can overshoot
+// the per-call limit by up to one step. consumed must be capped at the number
+// of fully available input positions or the history trim is silently skipped
+// while the accumulator is rebased, corrupting output and leaking memory.
+func TestPolyphase_SevereDownsampling_MonotonicAndBounded(t *testing.T) {
+	cases := []struct{ in, out float64 }{
+		{48000, 3001},
+		{48000, 1000.5},
+		{48000, 200.5},
+	}
+	for _, c := range cases {
+		r, err := NewResampler[float64](c.in, c.out, QualityHigh)
+		if err != nil {
+			t.Fatalf("%v to %v: %v", c.in, c.out, err)
+		}
+		const (
+			chunk = 4800
+			calls = 200
+			// transientSkip ignores the monotonic check over the initial
+			// filter transient, where priming legitimately dips before the
+			// delay line fills.
+			transientSkip = 100
+			// monotonicEps tolerates float rounding when comparing adjacent
+			// ramp outputs for non-decreasing order.
+			monotonicEps = 1e-6
+			// countSlackHigh/Low bound how far the total output may exceed or
+			// fall short of the ideal n*ratio: the high side absorbs
+			// boundary-carry rounding across many chunks, the low side absorbs
+			// the startup latency withheld at these severe ratios.
+			countSlackHigh = 64
+			countSlackLow  = 256
+			// historyChunkSlack bounds the retained delay line: steady-state
+			// history never exceeds tapsPerPhase-1 plus a few chunks of slack.
+			historyChunkSlack = chunk * 4
+		)
+		x := 0.0
+		last := -1.0
+		total := 0
+		for call := range calls {
+			in := make([]float64, chunk)
+			for i := range in {
+				in[i] = x
+				x += 1.0
+			}
+			out, err := r.Process(in)
+			if err != nil {
+				t.Fatalf("%v to %v call %d: %v", c.in, c.out, call, err)
+			}
+			total += len(out)
+			for i, v := range out {
+				// Ramp input must produce non-decreasing output away from
+				// the initial filter transient.
+				if total > transientSkip && v < last-monotonicEps {
+					t.Fatalf("%v to %v call %d sample %d: non-monotonic %g after %g",
+						c.in, c.out, call, i, v, last)
+				}
+				last = v
+			}
+		}
+		ratio := c.out / c.in
+		expected := float64(calls*chunk) * ratio
+		if float64(total) > expected+countSlackHigh || float64(total) < expected-countSlackLow {
+			t.Fatalf("%v to %v: total output %d, expected about %.0f", c.in, c.out, total, expected)
+		}
+		if r.polyphaseStage != nil {
+			maxHist := r.polyphaseStage.tapsPerPhase - 1 + historyChunkSlack
+			if len(r.polyphaseStage.history) > maxHist {
+				t.Fatalf("%v to %v: history grew to %d (bound %d)", c.in, c.out,
+					len(r.polyphaseStage.history), maxHist)
+			}
+		}
+	}
+}
diff --git a/internal/engine/stage_adapter.go b/internal/engine/stage_adapter.go
index a3adcd6..0ba9580 100644
--- a/internal/engine/stage_adapter.go
+++ b/internal/engine/stage_adapter.go
@@ -40,6 +40,11 @@ func (s *StageAdapter[F]) GetRatio() float64 {
 
 // GetLatency returns the stage latency in samples.
 // This is the delay due to FIR filter buffering.
+//
+// This is a filter group-delay heuristic in the input domain, consumed by
+// GetInfo for reporting. Streaming users priming a FIFO with silence should
+// instead use the engine Resampler.Latency() accessor, which reports the
+// startup deficit in output samples.
 func (s *StageAdapter[F]) GetLatency() int {
 	latency := 0
 
@@ -53,6 +58,16 @@ func (s *StageAdapter[F]) GetLatency() int {
 		latency += s.polyphaseStage.tapsPerPhase / latencyDivisor
 	}
 
+	// DFT decimation stage latency
+	if s.decimationStage != nil {
+		latency += s.decimationStage.numTaps / latencyDivisor
+	}
+
+	// Cubic interpolation stage latency
+	if s.cubicStage != nil {
+		latency += cubicLatencySamples
+	}
+
 	return latency
 }
 
@@ -84,13 +99,33 @@ func (s *StageAdapter[F]) GetMemoryUsage() int64 {
 
 	// Polyphase stage memory
 	if s.polyphaseStage != nil {
-		// Phase-first layout: polyCoeffs[phase][tap]
-		for _, phase := range s.polyphaseStage.polyCoeffs {
-			usage += int64(len(phase)) * bytesPerElement
+		// Phase-first layout: polyCoeffs[phase][tap]. All four cubic-interpolation
+		// coefficient banks (a, b, c, d) are the same shape and must all be
+		// counted; summing only polyCoeffs undercounts the coefficients by 4x.
+		for _, bank := range [][][]F{
+			s.polyphaseStage.polyCoeffs,
+			s.polyphaseStage.polyCoeffsB,
+			s.polyphaseStage.polyCoeffsC,
+			s.polyphaseStage.polyCoeffsD,
+		} {
+			for _, phase := range bank {
+				usage += int64(len(phase)) * bytesPerElement
+			}
 		}
 		usage += int64(cap(s.polyphaseStage.history)) * bytesPerElement
 	}
 
+	// Decimation stage memory
+	if s.decimationStage != nil {
+		usage += int64(len(s.decimationStage.coeffs)) * bytesPerElement
+		usage += int64(cap(s.decimationStage.history)) * bytesPerElement
+	}
+
+	// Cubic interpolation stage memory (symmetry with GetLatency).
+	if s.cubicStage != nil {
+		usage += s.cubicStage.GetMemoryUsage()
+	}
+
 	return usage
 }
 
diff --git a/internal/engine/stage_adapter_test.go b/internal/engine/stage_adapter_test.go
index e664f6a..160f7ea 100644
--- a/internal/engine/stage_adapter_test.go
+++ b/internal/engine/stage_adapter_test.go
@@ -32,3 +32,97 @@ func TestStageAdapterFloat32_ProcessZeroCopy(t *testing.T) {
 		t.Fatalf("ProcessZeroCopy returned error: %v", err)
 	}
 }
+
+// The public resample-package pipeline never wires QualityQuick through a
+// StageAdapter, so GetLatency's cubic branch is only reachable by constructing
+// the adapter directly. Pin it here.
+func TestStageAdapter_GetLatency_CubicBranch(t *testing.T) {
+	r, err := NewResampler[float64](44100, 48000, QualityQuick)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if r.cubicStage == nil {
+		t.Fatal("QualityQuick did not create a cubic stage")
+	}
+	adapter := NewStageAdapter(r)
+	if got := adapter.GetLatency(); got != cubicLatencySamples {
+		t.Errorf("GetLatency() = %d, want %d (cubic branch)", got, cubicLatencySamples)
+	}
+}
+
+// A QualityQuick resampler has only a cubic stage, so GetMemoryUsage must
+// report exactly the cubic stage's own accounting.
+func TestStageAdapter_GetMemoryUsage_CubicBranch(t *testing.T) {
+	r, err := NewResampler[float64](44100, 48000, QualityQuick)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if r.cubicStage == nil {
+		t.Fatal("QualityQuick did not create a cubic stage")
+	}
+	adapter := NewStageAdapter(r)
+	want := r.cubicStage.GetMemoryUsage()
+	if got := adapter.GetMemoryUsage(); got != want {
+		t.Errorf("GetMemoryUsage() = %d, want %d (cubic branch)", got, want)
+	}
+}
+
+// Integer 2:1 downsampling wires only the DFT decimation stage; exercise the
+// decimation branch of GetMemoryUsage.
+func TestStageAdapter_GetMemoryUsage_DecimationBranch(t *testing.T) {
+	r, err := NewResampler[float64](96000, 48000, QualityMedium)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if r.decimationStage == nil {
+		t.Fatal("expected a decimation stage for 96000->48000")
+	}
+	adapter := NewStageAdapter(r)
+	const bytesPerElement = int64(bytesPerFloat64)
+	want := int64(len(r.decimationStage.coeffs))*bytesPerElement +
+		int64(cap(r.decimationStage.history))*bytesPerElement
+	if got := adapter.GetMemoryUsage(); got != want {
+		t.Errorf("GetMemoryUsage() = %d, want %d (decimation branch)", got, want)
+	}
+}
+
+// GetMemoryUsage must count all four polyphase coefficient banks (a, b, c, d),
+// not just the base bank; counting only polyCoeffs undercounts the
+// coefficients by 4x.
+func TestStageAdapter_GetMemoryUsage_CountsAllFourPolyphaseBanks(t *testing.T) {
+	// Non-integer downsampling wires a DFT pre-stage plus a polyphase stage
+	// whose four cubic-interpolation coefficient banks all consume memory.
+	r, err := NewResampler[float64](48000, 44100, QualityMedium)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if r.polyphaseStage == nil {
+		t.Fatal("expected a polyphase stage for 48000->44100")
+	}
+	adapter := NewStageAdapter(r)
+
+	const bytesPerElement = int64(bytesPerFloat64)
+	bankBytes := func(bank [][]float64) int64 {
+		var n int64
+		for _, phase := range bank {
+			n += int64(len(phase)) * bytesPerElement
+		}
+		return n
+	}
+
+	ps := r.polyphaseStage
+	oneBank := bankBytes(ps.polyCoeffs)
+	if oneBank == 0 {
+		t.Fatal("polyphase coefficient bank is empty; test cannot distinguish the undercount")
+	}
+	want := oneBank + bankBytes(ps.polyCoeffsB) + bankBytes(ps.polyCoeffsC) + bankBytes(ps.polyCoeffsD)
+	want += int64(cap(ps.history)) * bytesPerElement
+	if r.preStage != nil {
+		want += bankBytes(r.preStage.polyCoeffs)
+		want += int64(cap(r.preStage.history)) * bytesPerElement
+	}
+
+	if got := adapter.GetMemoryUsage(); got != want {
+		t.Errorf("GetMemoryUsage() = %d, want %d (all four polyphase banks must be counted)", got, want)
+	}
+}
diff --git a/latency_test.go b/latency_test.go
new file mode 100644
index 0000000..2a6430f
--- /dev/null
+++ b/latency_test.go
@@ -0,0 +1,70 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import (
+	"math"
+	"testing"
+)
+
+// Latency() must predict the startup deficit: how many output samples the
+// first Process call withholds while the filter primes. Verified against
+// the measured deficit for a large first chunk.
+func TestLatency_MatchesMeasuredDeficit(t *testing.T) {
+	for _, c := range []struct{ in, out float64 }{
+		{44100, 48000},
+		{48000, 44100},
+		{48000, 16000},
+		{48000, 48000},
+	} {
+		for _, q := range []QualityPreset{QualityQuick, QualityLow, QualityMedium, QualityHigh} {
+			r, err := NewEngine(c.in, c.out, q)
+			if err != nil {
+				t.Fatal(err)
+			}
+			const n = 44100
+			in := make([]float64, n)
+			for i := range in {
+				in[i] = 0.5 * math.Sin(2*math.Pi*997*float64(i)/c.in)
+			}
+			out, err := r.Process(in)
+			if err != nil {
+				t.Fatal(err)
+			}
+			measured := int(float64(n)*c.out/c.in) - len(out)
+			got := r.Latency()
+			if got < measured-2 || got > measured+2 {
+				t.Errorf("%v to %v q=%v: Latency()=%d, measured deficit %d",
+					c.in, c.out, q, got, measured)
+			}
+		}
+	}
+}
+
+// Float32 mirror of one row of TestLatency_MatchesMeasuredDeficit: the float32
+// engine (NewEngineFloat32) must report the same startup deficit as the
+// float64 path.
+func TestLatencyFloat32_MatchesMeasuredDeficit(t *testing.T) {
+	const in, out = 44100.0, 48000.0
+	for _, q := range []QualityPreset{QualityQuick, QualityLow, QualityMedium, QualityHigh} {
+		r, err := NewEngineFloat32(in, out, q)
+		if err != nil {
+			t.Fatal(err)
+		}
+		const n = 44100
+		sig := make([]float32, n)
+		for i := range sig {
+			sig[i] = float32(0.5 * math.Sin(2*math.Pi*997*float64(i)/in))
+		}
+		outSamples, err := r.Process(sig)
+		if err != nil {
+			t.Fatal(err)
+		}
+		measured := int(float64(n)*out/in) - len(outSamples)
+		got := r.Latency()
+		if got < measured-2 || got > measured+2 {
+			t.Errorf("q=%v: Latency()=%d, measured deficit %d", q, got, measured)
+		}
+	}
+}
diff --git a/nan_validation_test.go b/nan_validation_test.go
new file mode 100644
index 0000000..1a74b84
--- /dev/null
+++ b/nan_validation_test.go
@@ -0,0 +1,52 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import (
+	"math"
+	"testing"
+)
+
+func TestNewEngine_RejectsNaNRates(t *testing.T) {
+	nan := math.NaN()
+	for _, c := range []struct{ in, out float64 }{
+		{nan, 48000},
+		{48000, nan},
+		{nan, nan},
+	} {
+		if _, err := NewEngine(c.in, c.out, QualityHigh); err == nil {
+			t.Errorf("NewEngine(%v, %v) accepted NaN", c.in, c.out)
+		}
+		if _, err := NewEngineFloat32(c.in, c.out, QualityHigh); err == nil {
+			t.Errorf("NewEngineFloat32(%v, %v) accepted NaN", c.in, c.out)
+		}
+	}
+}
+
+// Config.Validate (and thus the pipeline constructors New/NewMultiChannel/
+// NewStereo/NewSimple and the preset helpers) must reject NaN rates. NaN-blind
+// comparisons (c.InputRate <= 0, ratio < min) silently pass every branch, so a
+// NaN rate would otherwise build a pipeline that produces unresampled
+// passthrough garbage instead of returning an error.
+func TestNew_RejectsNaNRates(t *testing.T) {
+	nan := math.NaN()
+	for _, c := range []struct{ in, out float64 }{
+		{nan, 48000},
+		{48000, nan},
+		{nan, nan},
+	} {
+		cfg := &Config{
+			InputRate:  c.in,
+			OutputRate: c.out,
+			Channels:   1,
+			Quality:    QualitySpec{Preset: QualityHigh},
+		}
+		if _, err := New(cfg); err == nil {
+			t.Errorf("New(InputRate=%v, OutputRate=%v) accepted NaN", c.in, c.out)
+		}
+		if _, err := NewMultiChannel(c.in, c.out, 2, QualityHigh); err == nil {
+			t.Errorf("NewMultiChannel(%v, %v) accepted NaN", c.in, c.out)
+		}
+	}
+}
diff --git a/pipeline_builder.go b/pipeline_builder.go
index c9b3e98..01825b4 100644
--- a/pipeline_builder.go
+++ b/pipeline_builder.go
@@ -81,7 +81,7 @@ func createStage(spec StageSpec, config *Config) (Stage, error) {
 		return newCubicStage(spec.Ratio), nil
 
 	case pipeline.StageHalfBand:
-		return newHalfBandStage(spec.Ratio, spec.FilterLength, precision), nil
+		return newHalfBandStage(spec.Ratio, spec.FilterLength, precision)
 
 	case pipeline.StagePolyphase:
 		return newPolyphaseStage(
diff --git a/processinto_test.go b/processinto_test.go
index 8fe8a47..3a310b8 100644
--- a/processinto_test.go
+++ b/processinto_test.go
@@ -9,6 +9,8 @@ import (
 	"math"
 	"math/rand"
 	"testing"
+
+	"github.com/stretchr/testify/require"
 )
 
 type processIntoResampler interface {
@@ -139,10 +141,12 @@ func TestProcessInto_ZeroAllocs(t *testing.T) { //nolint:dupl // intentional par
 				t.Fatal(err)
 			}
 
+			var runErr error
 			allocs := testing.AllocsPerRun(100, func() {
 				r.Reset()
-				_, _ = r.ProcessInto(input, output)
+				_, runErr = r.ProcessInto(input, output)
 			})
+			require.NoError(t, runErr)
 
 			if allocs != 0 {
 				t.Fatalf("ProcessInto allocated %.0f times per call; expected 0", allocs)
@@ -151,6 +155,65 @@ func TestProcessInto_ZeroAllocs(t *testing.T) { //nolint:dupl // intentional par
 	}
 }
 
+// assertProcessIntoWarmZeroAllocs warms the resampler once, then asserts that
+// a Reset+ProcessInto cycle allocates nothing under testing.AllocsPerRun.
+// Resetting inside the measured closure keeps the internal buffers at their
+// grown capacity while re-exercising the priming path each iteration.
+func assertProcessIntoWarmZeroAllocs(t *testing.T, r *SimpleResampler, input, output []float64) {
+	t.Helper()
+
+	r.Reset()
+	if _, err := r.ProcessInto(input, output); err != nil {
+		t.Fatal(err)
+	}
+
+	var runErr error
+	allocs := testing.AllocsPerRun(100, func() {
+		r.Reset()
+		_, runErr = r.ProcessInto(input, output)
+	})
+	require.NoError(t, runErr)
+	if allocs != 0 {
+		t.Fatalf("ProcessInto allocated %.0f times per call; expected 0", allocs)
+	}
+}
+
+// TestProcessInto_ZeroAllocs_QualityQuick enforces the zero-allocation
+// invariant for the QualityQuick (cubic) path. doc.go advertises ProcessInto
+// as zero-allocation once warm; on this branch QualityQuick maps to the cubic
+// stage, so that stage must reuse a persistent output buffer like the FIR
+// stages do instead of allocating a fresh slice per call. TestProcessInto_ZeroAllocs
+// covers only the FIR qualities.
+func TestProcessInto_ZeroAllocs_QualityQuick(t *testing.T) {
+	cases := []struct {
+		name            string
+		inRate, outRate float64
+		durSeconds      int
+	}{
+		{"44100to48_quick", 44100, 48000, 3}, // upsample
+		{"48to16_quick", 48000, 16000, 3},    // integer downsample
+		{"48to44100_quick", 48000, 44100, 3}, // non-integer downsample
+	}
+
+	for _, tc := range cases {
+		t.Run(tc.name, func(t *testing.T) {
+			r, err := NewEngine(tc.inRate, tc.outRate, QualityQuick)
+			if err != nil {
+				t.Fatal(err)
+			}
+
+			inputLen := int(tc.inRate) * tc.durSeconds
+			input := make([]float64, inputLen)
+			for i := range input {
+				input[i] = float64(i) * 1e-5
+			}
+			output := make([]float64, r.EstimateOutput(inputLen))
+
+			assertProcessIntoWarmZeroAllocs(t, r, input, output)
+		})
+	}
+}
+
 // TestProcessInto_BufferTooSmall verifies that ErrBufferTooSmall is returned
 // when the output buffer is insufficient.
 func TestProcessInto_BufferTooSmall(t *testing.T) {
@@ -484,10 +547,12 @@ func TestProcessIntoFloat32_ZeroAllocs(t *testing.T) { //nolint:dupl // intentio
 				t.Fatal(err)
 			}
 
+			var runErr error
 			allocs := testing.AllocsPerRun(100, func() {
 				r.Reset()
-				_, _ = r.ProcessInto(input, output)
+				_, runErr = r.ProcessInto(input, output)
 			})
+			require.NoError(t, runErr)
 
 			if allocs != 0 {
 				t.Fatalf("ProcessInto allocated %.0f times per call; expected 0", allocs)
@@ -496,6 +561,63 @@ func TestProcessIntoFloat32_ZeroAllocs(t *testing.T) { //nolint:dupl // intentio
 	}
 }
 
+// assertProcessIntoFloat32WarmZeroAllocs is the float32 mirror of
+// assertProcessIntoWarmZeroAllocs: warm the engine once, then assert a
+// Reset+ProcessInto cycle allocates nothing under testing.AllocsPerRun.
+func assertProcessIntoFloat32WarmZeroAllocs(t *testing.T, r *SimpleResamplerFloat32, input, output []float32) {
+	t.Helper()
+
+	r.Reset()
+	if _, err := r.ProcessInto(input, output); err != nil {
+		t.Fatal(err)
+	}
+
+	var runErr error
+	allocs := testing.AllocsPerRun(100, func() {
+		r.Reset()
+		_, runErr = r.ProcessInto(input, output)
+	})
+	require.NoError(t, runErr)
+	if allocs != 0 {
+		t.Fatalf("ProcessInto allocated %.0f times per call; expected 0", allocs)
+	}
+}
+
+// TestProcessIntoFloat32_ZeroAllocs_QualityQuick enforces the zero-allocation
+// invariant for the QualityQuick (cubic) path on the float32 engine, mirroring
+// TestProcessInto_ZeroAllocs_QualityQuick. QualityQuick maps to the cubic
+// stage, which must reuse a persistent output buffer instead of allocating a
+// fresh slice per call.
+func TestProcessIntoFloat32_ZeroAllocs_QualityQuick(t *testing.T) {
+	cases := []struct {
+		name            string
+		inRate, outRate float64
+		durSeconds      int
+	}{
+		{"44100to48_quick", 44100, 48000, 3}, // upsample
+		{"48to16_quick", 48000, 16000, 3},    // integer downsample
+		{"48to44100_quick", 48000, 44100, 3}, // non-integer downsample
+	}
+
+	for _, tc := range cases {
+		t.Run(tc.name, func(t *testing.T) {
+			r, err := NewEngineFloat32(tc.inRate, tc.outRate, QualityQuick)
+			if err != nil {
+				t.Fatal(err)
+			}
+
+			inputLen := int(tc.inRate) * tc.durSeconds
+			input := make([]float32, inputLen)
+			for i := range input {
+				input[i] = float32(i) * 1e-5
+			}
+			output := make([]float32, r.EstimateOutput(inputLen))
+
+			assertProcessIntoFloat32WarmZeroAllocs(t, r, input, output)
+		})
+	}
+}
+
 // TestProcessFloat32Into_ZeroAllocs enforces the zero-allocation invariant for
 // the float32 New(...) batch path (constantRateResampler.ProcessFloat32Into).
 // This path converts float32<->float64 through grow-only scratch buffers, so it
@@ -542,9 +664,11 @@ func TestProcessFloat32Into_ZeroAllocs(t *testing.T) {
 				t.Fatal(err)
 			}
 
+			var runErr error
 			allocs := testing.AllocsPerRun(100, func() {
-				_, _ = r.ProcessFloat32Into(input, output)
+				_, runErr = r.ProcessFloat32Into(input, output)
 			})
+			require.NoError(t, runErr)
 
 			if allocs != 0 {
 				t.Fatalf("ProcessFloat32Into allocated %.0f times per call; expected 0", allocs)
diff --git a/quality_quick_test.go b/quality_quick_test.go
new file mode 100644
index 0000000..a5804cf
--- /dev/null
+++ b/quality_quick_test.go
@@ -0,0 +1,20 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import "testing"
+
+// QualityQuick means cubic interpolation (resample.go documents this and
+// New() honors it); the NewEngine path must agree instead of silently
+// substituting a full FIR pipeline with different latency and flush
+// behavior.
+func TestNewEngine_QualityQuick_IsCubic(t *testing.T) {
+	r, err := NewEngine(44100, 48000, QualityQuick)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if got := r.Latency(); got > 8 {
+		t.Errorf("QualityQuick latency %d, want cubic-scale (<= 8): engine mapped to FIR pipeline", got)
+	}
+}
diff --git a/resample.go b/resample.go
index d9c0558..c8198cc 100644
--- a/resample.go
+++ b/resample.go
@@ -166,7 +166,11 @@ var (
 
 // Validate checks if the configuration is valid.
 func (c *Config) Validate() error {
-	if c.InputRate <= 0 || c.OutputRate <= 0 {
+	// Expressed positively (!(x > 0)) so a NaN rate fails validation instead of
+	// silently passing the <= 0 comparison, which every NaN comparison returns
+	// false for. A NaN rate would otherwise build a pipeline that produces
+	// unresampled passthrough garbage.
+	if !(c.InputRate > 0) || !(c.OutputRate > 0) {
 		return fmt.Errorf("%w: sample rates must be positive", ErrInvalidConfig)
 	}
 
@@ -178,8 +182,11 @@ func (c *Config) Validate() error {
 		return fmt.Errorf("%w: too many channels (max %d)", ErrInvalidConfig, maxChannels)
 	}
 
+	// Expressed positively so a NaN ratio (e.g. from Inf/Inf) fails the bounds
+	// check instead of silently passing every comparison, matching the guard
+	// in internal/engine/resampler.go.
 	ratio := c.OutputRate / c.InputRate
-	if ratio < minRatioFactor || ratio > maxRatioFactor {
+	if !(ratio >= minRatioFactor && ratio <= maxRatioFactor) {
 		return fmt.Errorf("%w: resampling ratio out of range (%v to %v)", ErrInvalidConfig, minRatioFactor, maxRatioFactor)
 	}
 
diff --git a/stages.go b/stages.go
index dedb721..69cfd77 100644
--- a/stages.go
+++ b/stages.go
@@ -4,6 +4,8 @@
 package resampler
 
 import (
+	"fmt"
+
 	"github.com/tphakala/go-audio-resampler/internal/engine"
 	"github.com/tphakala/go-audio-resampler/internal/pipeline"
 )
@@ -28,19 +30,22 @@ func newCubicStage(ratio float64) pipeline.Stage {
 // For now, we use the polyphase engine which is functionally equivalent but less optimized.
 // A dedicated half-band implementation would store only non-zero coefficients
 // (see soxr/src/half-coefs.h for reference coefficients).
-func newHalfBandStage(ratio float64, filterLength, precision int) pipeline.Stage {
+//
+// The underlying polyphase construction can fail in principle (invalid filter
+// parameters), so the error is propagated rather than silently substituted
+// with a degraded stub. pipeline.BuildPipeline only ever emits StageHalfBand
+// specs with Ratio fixed at halfRatio (0.5) or doubleRatio (2.0), both well
+// inside engine.NewResampler's valid ratio range, so construction failure is
+// not reachable through any public New() input today; propagating the error
+// keeps that guarantee explicit instead of relying on silent substitution.
+func newHalfBandStage(ratio float64, filterLength, precision int) (pipeline.Stage, error) {
 	// Use polyphase stage with factor=2 - functionally equivalent to half-band
 	// The polyphase implementation handles 2x ratios efficiently
 	stage, err := newPolyphaseStage(ratio, filterLength, halfBandFactor, precision)
 	if err != nil {
-		// Fallback to stub if polyphase creation fails
-		return &stubStage{
-			ratio:        ratio,
-			filterLength: filterLength,
-			name:         "halfband",
-		}
+		return nil, fmt.Errorf("failed to create half-band stage: %w", err)
 	}
-	return stage
+	return stage, nil
 }
 
 // newPolyphaseStage creates a polyphase FIR filtering stage using engine.Resampler.
@@ -118,7 +123,9 @@ func newFFTStage(ratio float64, fftSize, precision int) (pipeline.Stage, error)
 	return newPolyphaseStage(ratio, fftSize, defaultFFTPhases, precision)
 }
 
-// stubStage is a temporary stub implementation for stages not yet implemented.
+// stubStage is a minimal passthrough Stage used only by tests; no production
+// code path constructs it. It nearest-neighbor resamples by the ratio, which
+// is enough to exercise the pipeline.Stage interface without a real filter.
 type stubStage struct {
 	ratio        float64
 	filterLength int
diff --git a/stages_test.go b/stages_test.go
new file mode 100644
index 0000000..99a2583
--- /dev/null
+++ b/stages_test.go
@@ -0,0 +1,139 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import (
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+	"github.com/stretchr/testify/require"
+
+	"github.com/tphakala/go-audio-resampler/internal/pipeline"
+)
+
+// newHalfBandStage must never substitute the degraded stubStage fallback for
+// a construction failure; it must return the error instead. This was
+// previously silent: a failed newPolyphaseStage call fell back to a
+// nearest-neighbor stub with no signal to the caller.
+//
+// pipeline.BuildPipeline only ever emits StageHalfBand specs with Ratio
+// pinned to halfRatio (0.5) or doubleRatio (2.0); those are the only two
+// ratios newHalfBandStage is ever invoked with through the public New()
+// path. Both are well inside engine.NewResampler's valid ratio range
+// [1/256, 256], so construction cannot fail for them today. These tests
+// pin the signature contract (error-free construction still works, and the
+// error return exists and would propagate if it were ever non-nil).
+func TestNewHalfBandStage_RepresentativeRatios(t *testing.T) {
+	ratios := []float64{0.5, 2.0} // exactly what pipeline.go emits for StageHalfBand
+	precisions := []int{16, 20, 24, 28, 32, 33}
+
+	for _, ratio := range ratios {
+		for _, precision := range precisions {
+			stage, err := newHalfBandStage(ratio, 0, precision)
+			require.NoErrorf(t, err, "ratio=%v precision=%d", ratio, precision)
+			require.NotNilf(t, stage, "ratio=%v precision=%d", ratio, precision)
+
+			// Must not be the degraded nearest-neighbor stub.
+			if _, isStub := stage.(*stubStage); isStub {
+				t.Fatalf("ratio=%v precision=%d: newHalfBandStage returned stubStage instead of a real polyphase stage", ratio, precision)
+			}
+
+			assert.InDeltaf(t, ratio, stage.GetRatio(), 1e-9, "ratio=%v precision=%d", ratio, precision)
+		}
+	}
+}
+
+// createStage's StageHalfBand branch must propagate newHalfBandStage's
+// (Stage, error) return directly, matching the pattern already used by the
+// StagePolyphase and StageFFT branches, instead of discarding the error and
+// always returning nil.
+func TestCreateStage_HalfBand_PropagatesSignature(t *testing.T) {
+	config := &Config{
+		InputRate:  48000,
+		OutputRate: 24000,
+		Channels:   1,
+		Quality:    GetPresetSpec(QualityHigh),
+	}
+
+	spec := StageSpec{
+		StageSpec: pipeline.StageSpec{
+			Type:         pipeline.StageHalfBand,
+			Ratio:        0.5,
+			FilterLength: 32,
+		},
+	}
+
+	stage, err := createStage(spec, config)
+	require.NoError(t, err)
+	require.NotNil(t, stage)
+
+	if _, isStub := stage.(*stubStage); isStub {
+		t.Fatal("createStage(StageHalfBand) returned stubStage instead of a real polyphase stage")
+	}
+	assert.InDelta(t, 0.5, stage.GetRatio(), 1e-9)
+}
+
+// New() must construct a working pipeline (no silently degraded stub
+// stages) for configurations whose ratio decomposition routes through one
+// or more half-band stages, in both the upsampling and downsampling
+// directions.
+func TestNew_HalfBandPipeline_NoStubFallback(t *testing.T) {
+	cases := []struct {
+		name       string
+		inputRate  float64
+		outputRate float64
+	}{
+		{"upsample_needs_halfband", 8000, 48000},   // ratio 6.0: halfband stages factor out powers of 2
+		{"downsample_needs_halfband", 48000, 8000}, // ratio 1/6: halfband stages factor out powers of 2
+	}
+
+	for _, tc := range cases {
+		t.Run(tc.name, func(t *testing.T) {
+			config := &Config{
+				InputRate:  tc.inputRate,
+				OutputRate: tc.outputRate,
+				Channels:   1,
+				Quality:    QualitySpec{Preset: QualityHigh},
+			}
+
+			r, err := New(config)
+			require.NoError(t, err)
+			require.NotNil(t, r)
+
+			crr, ok := r.(*constantRateResampler)
+			require.True(t, ok, "expected *constantRateResampler")
+
+			sawHalfBand := false
+			for _, spec := range crr.pipeline.stages {
+				if spec.engine == "halfband" {
+					sawHalfBand = true
+				}
+			}
+			require.True(t, sawHalfBand, "test config expected to produce at least one halfband stage")
+
+			for chIdx, ch := range crr.channels {
+				for stIdx, stg := range ch.stages {
+					if _, isStub := stg.(*stubStage); isStub {
+						t.Fatalf("channel %d stage %d is a stubStage: half-band construction silently degraded", chIdx, stIdx)
+					}
+				}
+			}
+
+			// Confirm the pipeline actually processes audio end to end.
+			// Input must clear filter startup latency to produce output from
+			// Process alone, so also drain Flush to be independent of that.
+			input := make([]float64, 8192)
+			for i := range input {
+				input[i] = 1.0
+			}
+			out, err := r.Process(input)
+			require.NoError(t, err)
+
+			tail, err := r.Flush()
+			require.NoError(t, err)
+
+			assert.NotEmpty(t, append(out, tail...))
+		})
+	}
+}
diff --git a/streaming_equivalence_test.go b/streaming_equivalence_test.go
new file mode 100644
index 0000000..93b5e6b
--- /dev/null
+++ b/streaming_equivalence_test.go
@@ -0,0 +1,181 @@
+// SPDX-FileCopyrightText: 2026 Tomi P. Hakala
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+package resampler
+
+import (
+	"math"
+	"math/rand"
+	"slices"
+	"testing"
+)
+
+// Chunked streaming Process calls plus one final Flush must be bit-exact
+// with processing the whole signal in a single Process+Flush. This is the
+// contract that makes real-time chunked use (issue #51) safe.
+func TestStreamingEquivalence_Float64(t *testing.T) {
+	ratios := []struct {
+		name    string
+		in, out float64
+	}{
+		{"44k1_to_48k", 44100, 48000},
+		{"48k_to_44k1", 48000, 44100},
+		{"48k_to_16k_integer", 48000, 16000},
+		{"unity", 48000, 48000},
+	}
+	qualities := []QualityPreset{QualityLow, QualityMedium, QualityHigh}
+	chunkPlans := [][]int{
+		{1, 7, 13, 470, 4096},
+		{31, 331, 997},
+	}
+
+	const n = 44100
+	input := make([]float64, n)
+	for i := range input {
+		input[i] = 0.5 * math.Sin(2*math.Pi*997*float64(i)/44100)
+	}
+
+	for _, rr := range ratios {
+		for _, q := range qualities {
+			oneShot, err := NewEngine(rr.in, rr.out, q)
+			if err != nil {
+				t.Fatalf("%s: NewEngine: %v", rr.name, err)
+			}
+			ref, err := oneShot.Process(slices.Clone(input))
+			if err != nil {
+				t.Fatalf("%s: Process: %v", rr.name, err)
+			}
+			refTail, err := oneShot.Flush()
+			if err != nil {
+				t.Fatalf("%s: Flush: %v", rr.name, err)
+			}
+			ref = append(ref, refTail...)
+
+			for pi, plan := range chunkPlans {
+				assertChunkedPlanEqualsFloat64(t, rr.name, rr.in, rr.out, q, pi, plan, input, ref)
+			}
+		}
+	}
+}
+
+// assertChunkedPlanEqualsFloat64 feeds input through a fresh engine in
+// pseudo-random chunk sizes drawn from plan (seeded by pi for reproducibility),
+// Flushes once, and asserts the result is bit-exact with ref.
+func assertChunkedPlanEqualsFloat64(t *testing.T, name string, in, out float64, q QualityPreset, pi int, plan []int, input, ref []float64) {
+	t.Helper()
+
+	chunked, err := NewEngine(in, out, q)
+	if err != nil {
+		t.Fatalf("%s: NewEngine: %v", name, err)
+	}
+	var got []float64
+	rng := rand.New(rand.NewSource(int64(pi) + 1))
+	for pos := 0; pos < len(input); {
+		size := plan[rng.Intn(len(plan))]
+		if pos+size > len(input) {
+			size = len(input) - pos
+		}
+		outChunk, err := chunked.Process(slices.Clone(input[pos : pos+size]))
+		if err != nil {
+			t.Fatalf("%s plan %d: Process: %v", name, pi, err)
+		}
+		got = append(got, outChunk...)
+		pos += size
+	}
+	tail, err := chunked.Flush()
+	if err != nil {
+		t.Fatalf("%s plan %d: Flush: %v", name, pi, err)
+	}
+	got = append(got, tail...)
+
+	if len(got) != len(ref) {
+		t.Fatalf("%s q=%v plan %d: length %d != one-shot %d", name, q, pi, len(got), len(ref))
+	}
+	for i := range got {
+		if got[i] != ref[i] {
+			t.Fatalf("%s q=%v plan %d: sample %d differs: %g != %g", name, q, pi, i, got[i], ref[i])
+		}
+	}
+}
+
+func TestStreamingEquivalence_Float32(t *testing.T) {
+	// Broadened beyond the single issue #51 configuration: an upsample and a
+	// downsample ratio, two qualities, and two fixed chunk sizes. Chunked
+	// Process plus one Flush must stay bit-exact with the one-shot path.
+	const n = 44100
+	ratios := []struct {
+		name    string
+		in, out float64
+	}{
+		{"44k1_to_48k", 44100, 48000},
+		{"48k_to_44k1", 48000, 44100},
+	}
+	qualities := []QualityPreset{QualityMedium, QualityHigh}
+	chunkSizes := []int{470, 997}
+
+	for _, rr := range ratios {
+		input := make([]float32, n)
+		for i := range input {
+			input[i] = float32(0.5 * math.Sin(2*math.Pi*997*float64(i)/rr.in))
+		}
+		for _, q := range qualities {
+			for _, size := range chunkSizes {
+				assertStreamingEquivalentFloat32(t, rr.name, rr.in, rr.out, q, size, input)
+			}
+		}
+	}
+}
+
+// assertStreamingEquivalentFloat32 checks that feeding input in fixed-size
+// chunks (with one final Flush) is bit-exact with a single one-shot
+// Process+Flush for the float32 engine at the given ratio, quality, and chunk
+// size.
+func assertStreamingEquivalentFloat32(t *testing.T, name string, in, out float64, q QualityPreset, size int, input []float32) {
+	t.Helper()
+
+	oneShot, err := NewEngineFloat32(in, out, q)
+	if err != nil {
+		t.Fatalf("%s q=%v: NewEngineFloat32: %v", name, q, err)
+	}
+	ref, err := oneShot.Process(slices.Clone(input))
+	if err != nil {
+		t.Fatalf("%s q=%v: Process: %v", name, q, err)
+	}
+	refTail, err := oneShot.Flush()
+	if err != nil {
+		t.Fatalf("%s q=%v: Flush: %v", name, q, err)
+	}
+	ref = append(ref, refTail...)
+
+	chunked, err := NewEngineFloat32(in, out, q)
+	if err != nil {
+		t.Fatalf("%s q=%v size=%d: NewEngineFloat32: %v", name, q, size, err)
+	}
+	var got []float32
+	for pos := 0; pos < len(input); {
+		step := size
+		if pos+step > len(input) {
+			step = len(input) - pos
+		}
+		outChunk, err := chunked.Process(slices.Clone(input[pos : pos+step]))
+		if err != nil {
+			t.Fatalf("%s q=%v size=%d: Process: %v", name, q, size, err)
+		}
+		got = append(got, outChunk...)
+		pos += step
+	}
+	tail, err := chunked.Flush()
+	if err != nil {
+		t.Fatalf("%s q=%v size=%d: Flush: %v", name, q, size, err)
+	}
+	got = append(got, tail...)
+
+	if len(got) != len(ref) {
+		t.Fatalf("%s q=%v size=%d: length %d != one-shot %d", name, q, size, len(got), len(ref))
+	}
+	for i := range got {
+		if got[i] != ref[i] {
+			t.Fatalf("%s q=%v size=%d: sample %d differs: %g != %g", name, q, size, i, got[i], ref[i])
+		}
+	}
+}