From 6f983dbb055b3dba2e0e2fe18b79453e45919785 Mon Sep 17 00:00:00 2001
From: "Tomi P. Hakala" <7030001+tphakala@users.noreply.github.com>
Date: Fri, 17 Jul 2026 20:00:13 +0300
Subject: [PATCH 1/6] refactor(engine): extract shared cubic interp bank
 builder

The phase-wrap measurement test previously re-derived the Catmull-Rom
bank construction to build its differential reference, so a formula
change would have required lockstep edits. NewPolyphaseStage and the
test now share buildCubicInterpBanks and the test varies only the
boundary-indexing policy. The bank-equality assertion pins the
extraction bit-for-bit.

Refs #53
---
 internal/engine/phase_wrap_measure_test.go | 47 +++----------
 internal/engine/polyphase_stage.go         | 79 ++++++++++++----------
 2 files changed, 55 insertions(+), 71 deletions(-)

diff --git a/internal/engine/phase_wrap_measure_test.go b/internal/engine/phase_wrap_measure_test.go
index 06a46c8..df30903 100644
--- a/internal/engine/phase_wrap_measure_test.go
+++ b/internal/engine/phase_wrap_measure_test.go
@@ -78,14 +78,16 @@ type phaseWrapResult struct {
 	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.
+// buildPhaseInterpBanks builds the reference interpolation banks through the
+// shared production builder (buildCubicInterpBanks in polyphase_stage.go),
+// varying only the boundary-indexing policy. With flat=true it reproduces the
+// production (fixed) indexing: the neighbour of a phase boundary is the
+// adjacent prototype sample, crossing into the next/previous tap. With
+// flat=false it reproduces the original wrap: the neighbour phase is kept
+// within the same tap (phase % numPhases). Everything else, the Catmull-Rom
+// cubic math and the reversed tap storage, is identical (shared code), 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 {
@@ -109,36 +111,7 @@ func buildPhaseInterpBanks(fb *polyphaseFilter, numPhases int, flat bool) (a, b,
 		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
+	return buildCubicInterpBanks[float64](numPhases, fb.tapsPerPhase, getCoeff)
 }
 
 // processAllPhaseWrap runs Process followed by the terminal Flush and returns
diff --git a/internal/engine/polyphase_stage.go b/internal/engine/polyphase_stage.go
index 497d0e9..5c8d33b 100644
--- a/internal/engine/polyphase_stage.go
+++ b/internal/engine/polyphase_stage.go
@@ -57,6 +57,50 @@ type PolyphaseStage[F simdops.Float] struct {
 	samplesOut int64
 }
 
+// buildCubicInterpBanks builds the four Catmull-Rom sub-phase interpolation
+// coefficient banks from a prototype coefficient accessor. getCoeff(phase, tap)
+// returns the prototype coefficient for that phase and tap and applies the
+// caller's boundary-indexing policy for out-of-range phases. Taps are stored
+// reversed for the dot-product kernels. Shared by NewPolyphaseStage and the
+// phase-wrap measurement test so the interpolation formula lives in one place.
+func buildCubicInterpBanks[F simdops.Float](numPhases, tapsPerPhase int, getCoeff func(phase, tap int) float64) (a, b, c, d [][]F) {
+	a = make([][]F, numPhases)
+	b = make([][]F, numPhases)
+	c = make([][]F, numPhases)
+	d = make([][]F, numPhases)
+
+	for phase := range numPhases {
+		a[phase] = make([]F, tapsPerPhase)
+		b[phase] = make([]F, tapsPerPhase)
+		c[phase] = make([]F, tapsPerPhase)
+		d[phase] = make([]F, tapsPerPhase)
+
+		for tap := range tapsPerPhase {
+			// Get coefficients from adjacent phases for cubic interpolation
+			// f0 = current phase, f1 = next phase, fm1 = previous phase, f2 = next-next phase
+			f0 := getCoeff(phase, tap)
+			f1 := getCoeff(phase+1, tap)
+			fm1 := getCoeff(phase-1, tap)
+			f2 := getCoeff(phase+cubicPhaseOffset, tap)
+
+			// Compute cubic interpolation coefficients (Catmull-Rom style)
+			// These allow smooth interpolation: f(x) = a + b*x + c*x² + d*x³
+			av := f0
+			cv := cubicCenterCoeff*(f1+fm1) - f0
+			dv := (1.0 / cubicDivisor) * (f2 - f1 + fm1 - f0 - cubicCMultiplier*cv)
+			bv := f1 - f0 - dv - cv
+
+			// Store in REVERSED order for correct convolution direction
+			revTap := tapsPerPhase - 1 - tap
+			a[phase][revTap] = F(av)
+			b[phase][revTap] = F(bv)
+			c[phase][revTap] = F(cv)
+			d[phase][revTap] = F(dv)
+		}
+	}
+	return a, b, c, d
+}
+
 // NewPolyphaseStage creates a polyphase resampling stage.
 //
 // Parameters:
@@ -132,40 +176,7 @@ func NewPolyphaseStage[F simdops.Float](ratio, totalIORatio float64, hasPreStage
 	// Allocate coefficient arrays with cubic interpolation support
 	// polyCoeffs = a (base), polyCoeffsB = b (linear), polyCoeffsC = c (quadratic), polyCoeffsD = d (cubic)
 	// Interpolation formula: coef(x) = a + x*(b + x*(c + x*d)) where x ∈ [0, 1)
-	polyCoeffs := make([][]F, numPhases)
-	polyCoeffsB := make([][]F, numPhases)
-	polyCoeffsC := make([][]F, numPhases)
-	polyCoeffsD := make([][]F, numPhases)
-
-	for phase := range numPhases {
-		polyCoeffs[phase] = make([]F, tapsPerPhase)
-		polyCoeffsB[phase] = make([]F, tapsPerPhase)
-		polyCoeffsC[phase] = make([]F, tapsPerPhase)
-		polyCoeffsD[phase] = make([]F, tapsPerPhase)
-
-		for tap := range tapsPerPhase {
-			// Get coefficients from adjacent phases for cubic interpolation
-			// f0 = current phase, f1 = next phase, fm1 = previous phase, f2 = next-next phase
-			f0 := getCoeff(phase, tap)
-			f1 := getCoeff(phase+1, tap)
-			fm1 := getCoeff(phase-1, tap)
-			f2 := getCoeff(phase+cubicPhaseOffset, tap)
-
-			// Compute cubic interpolation coefficients (Catmull-Rom style)
-			// These allow smooth interpolation: f(x) = a + b*x + c*x² + d*x³
-			a := f0
-			c := cubicCenterCoeff*(f1+fm1) - f0
-			d := (1.0 / cubicDivisor) * (f2 - f1 + fm1 - f0 - cubicCMultiplier*c)
-			b := f1 - f0 - d - c
-
-			// Store in REVERSED order for correct convolution direction
-			revTap := tapsPerPhase - 1 - tap
-			polyCoeffs[phase][revTap] = F(a)
-			polyCoeffsB[phase][revTap] = F(b)
-			polyCoeffsC[phase][revTap] = F(c)
-			polyCoeffsD[phase][revTap] = F(d)
-		}
-	}
+	polyCoeffs, polyCoeffsB, polyCoeffsC, polyCoeffsD := buildCubicInterpBanks[F](numPhases, tapsPerPhase, getCoeff)
 
 	return &PolyphaseStage[F]{
 		polyCoeffs:    polyCoeffs,

From 701944f438c8af626dd0b6c387c0c7a8c060a816 Mon Sep 17 00:00:00 2001
From: "Tomi P. Hakala" <7030001+tphakala@users.noreply.github.com>
Date: Fri, 17 Jul 2026 20:02:04 +0300
Subject: [PATCH 2/6] test(engine): pin full-pipeline THD at
 active-interpolation ratios

quality_regression_test.go previously exercised only exact-rational or
integer ratios where the sub-phase interpolation banks are never
consulted, so the phase-boundary indexing fix from #56 was guarded only
at stage level. Add THD assertions through the whole resampler at
32000->44100 and 44100->64000 for High, Medium, and Low. High gets a
dedicated tighter limit (-145 dB, worst measured -151.63 dB); Medium and
Low fall within margin of the existing constants and reuse them.

Refs #54
---
 internal/engine/quality_regression_test.go | 47 +++++++++++++++++++++-
 1 file changed, 46 insertions(+), 1 deletion(-)

diff --git a/internal/engine/quality_regression_test.go b/internal/engine/quality_regression_test.go
index 1eda95b..54d36bc 100644
--- a/internal/engine/quality_regression_test.go
+++ b/internal/engine/quality_regression_test.go
@@ -44,6 +44,15 @@ const (
 	regressionMaxTHD_High     = -140.0 // dB (actual: ~-157 dB)
 	regressionMaxTHD_VeryHigh = -140.0 // dB (actual: ~-162 dB)
 
+	// Active-interpolation ratio THD limits (#54). Measured on 2026-07-17:
+	// 32000->44100 High: -162.05 dB, Medium: -134.95 dB, Low: -144.22 dB;
+	// 44100->64000 High: -151.63 dB, Medium: -135.03 dB, Low: -139.47 dB.
+	// Limits are measured values rounded up with ~5 dB margin. Medium and Low
+	// reuse the regressionMaxTHD_Medium and regressionMaxTHD_Low constants
+	// above since the measured values here fall within margin of those
+	// existing constants; only High needed a dedicated, tighter limit.
+	regressionMaxTHDInterp_High = -145.0 // dB (worst measured: -151.63 dB)
+
 	// Minimum SNR at 1kHz test frequency (more positive = better)
 	// Calibrated: actual SNR varies by conversion type (40-104 dB)
 	regressionMinSNR_Quick    = 35.0 // dB (actual: ~43 dB)
@@ -153,6 +162,43 @@ func TestQualityRegression_THD(t *testing.T) {
 	}
 }
 
+// Active-interpolation THD regression (#54): unlike the exact-rational
+// ratios above, these ratios consult the sub-phase Catmull-Rom
+// interpolation banks, so this pins the whole pipeline against the
+// phase-boundary indexing defect fixed in #56 (stage-level guard:
+// phase_wrap_measure_test.go). Thresholds are measured values plus margin;
+// see constants.
+func TestQualityRegression_THD_ActiveInterpolation(t *testing.T) {
+	tests := []struct {
+		inputRate, outputRate float64
+		testFreq              float64
+		quality               Quality
+		maxTHD                float64
+	}{
+		{32000, 44100, 1000, QualityHigh, regressionMaxTHDInterp_High},
+		{32000, 44100, 1000, QualityMedium, regressionMaxTHD_Medium},
+		{32000, 44100, 1000, QualityLow, regressionMaxTHD_Low},
+		{44100, 64000, 1000, QualityHigh, regressionMaxTHDInterp_High},
+		{44100, 64000, 1000, QualityMedium, regressionMaxTHD_Medium},
+		{44100, 64000, 1000, QualityLow, regressionMaxTHD_Low},
+	}
+
+	for _, tc := range tests {
+		name := qualityName(tc.quality) + "/" + formatRatio(tc.inputRate, tc.outputRate)
+		t.Run(name, func(t *testing.T) {
+			thd := measureTHDInternal(t, tc.inputRate, tc.outputRate, tc.testFreq, tc.quality)
+
+			// THD in dB: more negative is better
+			// If thd > maxTHD, it's worse (regression)
+			if thd > tc.maxTHD {
+				t.Errorf("THD REGRESSION: got %.2f dB, want <= %.2f dB", thd, tc.maxTHD)
+			} else {
+				t.Logf("THD: %.2f dB (threshold: %.2f dB) ✓", thd, tc.maxTHD)
+			}
+		})
+	}
+}
+
 // TestQualityRegression_SNR verifies SNR hasn't regressed
 func TestQualityRegression_SNR(t *testing.T) {
 	testCases := []struct {
@@ -501,4 +547,3 @@ func measurePassbandRippleInternal(t *testing.T, inputRate, outputRate float64,
 
 	return maxDev - minDev
 }
-

From 77e9c6011691ba5dae1e1a87f3ec53e956b4710d Mon Sep 17 00:00:00 2001
From: "Tomi P. Hakala" <7030001+tphakala@users.noreply.github.com>
Date: Fri, 17 Jul 2026 20:05:30 +0300
Subject: [PATCH 3/6] fix: report the true startup deficit from
 constantRateResampler.GetLatency

GetLatency summed per-stage group-delay heuristics whose terms live in
different rate domains, so multi-stage pipelines mis-reported the
startup deficit (672 reported vs 703 measured at 44100 to 96000
QualityHigh). The engine Resampler now exposes StartupDeficit(), the
un-rounded deficit in output samples that Latency() ceils, CubicStage
mirrors it, and GetLatency accumulates each stage's fractional deficit
through the downstream stage ratios before rounding once. Accuracy is
now within 2 samples across presets and ratios, pinned by a numeric
test on the New(config) path.

Refs #52
---
 CHANGELOG.md                 |  4 ++++
 constant.go                  | 31 ++++++++++++++++---------
 internal/engine/cubic.go     |  6 +++++
 internal/engine/resampler.go | 23 ++++++++++++------
 latency_test.go              | 45 ++++++++++++++++++++++++++++++++++++
 resample.go                  |  6 +++--
 6 files changed, 95 insertions(+), 20 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4752b68..58391a6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -43,6 +43,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - 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)
+- `GetLatency` on the `New(config)` path now reports the measured startup
+  deficit in output samples; it previously mixed rate domains across pipeline
+  stages and mis-reported multi-stage ratios (672 reported vs 703 measured at
+  44100 to 96000 QualityHigh; now within 2 samples). (#52)
 
 ### Changed
 
diff --git a/constant.go b/constant.go
index e6f7071..e8865fa 100644
--- a/constant.go
+++ b/constant.go
@@ -5,6 +5,7 @@ package resampler
 
 import (
 	"fmt"
+	"math"
 	"sync"
 
 	pipelinepkg "github.com/tphakala/go-audio-resampler/internal/pipeline"
@@ -403,26 +404,34 @@ func (r *constantRateResampler) FlushMulti() ([][]float64, error) {
 	return output, nil
 }
 
-// GetLatency returns the total pipeline latency in samples.
+// GetLatency returns the pipeline's startup deficit in output samples: how
+// many samples early Process calls withhold while the stage filters prime.
+// Each stage's deficit is converted into the final output rate domain
+// through the downstream stages' ratios before summing.
 func (r *constantRateResampler) GetLatency() int {
 	if r.pipeline == nil || len(r.channels) == 0 {
 		return 0
 	}
-
-	// Use the first channel's stages to calculate latency
 	ch := r.channels[0]
 	if ch == nil || len(ch.stages) == 0 {
 		return 0
 	}
-
-	totalLatency := 0
-	for _, stage := range ch.stages {
-		// Account for stage processing latency and ratio change
-		stageLatency := int(float64(stage.GetLatency()) * stage.GetRatio())
-		totalLatency += stageLatency
+	total := 0.0
+	for i, stage := range ch.stages {
+		var deficit float64
+		if s, ok := stage.(interface{ StartupDeficit() float64 }); ok {
+			deficit = s.StartupDeficit()
+		} else {
+			// Fallback for stages without deficit accounting: group-delay
+			// heuristic converted to the stage's output domain.
+			deficit = float64(stage.GetLatency()) * stage.GetRatio()
+		}
+		for _, downstream := range ch.stages[i+1:] {
+			deficit *= downstream.GetRatio()
+		}
+		total += deficit
 	}
-
-	return totalLatency
+	return int(math.Ceil(total))
 }
 
 // Reset clears all internal state.
diff --git a/internal/engine/cubic.go b/internal/engine/cubic.go
index ccc2830..0704fc7 100644
--- a/internal/engine/cubic.go
+++ b/internal/engine/cubic.go
@@ -192,6 +192,12 @@ func (c *CubicStage[F]) GetMemoryUsage() int64 {
 	return cubicMemoryUsage
 }
 
+// StartupDeficit returns the startup deficit in output samples, un-rounded.
+// Mirrors the cubic branch of Resampler.Latency for direct pipeline use.
+func (c *CubicStage[F]) StartupDeficit() float64 {
+	return cubicLatencySamples * c.ratio
+}
+
 // GetFilterLength returns 0 as cubic doesn't use a filter.
 func (c *CubicStage[F]) GetFilterLength() int {
 	return cubicInterpolationPoints
diff --git a/internal/engine/resampler.go b/internal/engine/resampler.go
index 75b9042..5c36d94 100644
--- a/internal/engine/resampler.go
+++ b/internal/engine/resampler.go
@@ -365,13 +365,14 @@ 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 {
+// StartupDeficit returns the startup deficit in output samples, un-rounded.
+// It exists so multi-stage pipelines can accumulate each stage's fractional
+// deficit and round once at the end, avoiding the per-stage rounding error
+// that accrues when every stage rounds to an integer before summing. Latency
+// wraps this with math.Ceil.
+func (r *Resampler[F]) StartupDeficit() float64 {
 	if r.cubicStage != nil {
-		return int(math.Ceil(cubicLatencySamples * r.ratio))
+		return cubicLatencySamples * r.ratio
 	}
 	deficitIn := 0.0
 	if r.preStage != nil && r.preStage.factor > 1 {
@@ -387,7 +388,15 @@ func (r *Resampler[F]) Latency() int {
 		}
 		deficitIn += float64(r.polyphaseStage.tapsPerPhase-1) / intermediateFactor
 	}
-	return int(math.Ceil(deficitIn * r.ratio))
+	return deficitIn * r.ratio
+}
+
+// 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 {
+	return int(math.Ceil(r.StartupDeficit()))
 }
 
 // isIntegerRatio checks if the ratio is an integer (within tolerance).
diff --git a/latency_test.go b/latency_test.go
index 2a6430f..9782c0d 100644
--- a/latency_test.go
+++ b/latency_test.go
@@ -68,3 +68,48 @@ func TestLatencyFloat32_MatchesMeasuredDeficit(t *testing.T) {
 		}
 	}
 }
+
+// GetLatency on the New(config) multi-stage path must report the startup
+// deficit in output samples, like SimpleResampler.Latency does for the
+// engine path. Guards against the domain-mixing error fixed in #52.
+func TestGetLatency_ConfigPath_MatchesMeasuredDeficit(t *testing.T) {
+	for _, c := range []struct{ in, out float64 }{
+		{44100, 48000},
+		{44100, 96000},
+		{48000, 44100},
+	} {
+		for _, preset := range []QualityPreset{QualityLow, QualityMedium, QualityHigh} {
+			cfg := Config{
+				InputRate:  c.in,
+				OutputRate: c.out,
+				Channels:   1,
+				Quality:    QualitySpec{Preset: preset},
+			}
+			r, err := New(&cfg)
+			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.GetLatency()
+			// The accumulated fractional deficit is rounded up once with
+			// math.Ceil, so GetLatency lands one or two samples above the
+			// measured deficit (observed max 2 across these rows). tol stays
+			// well below the old domain-mixing error, which under-reported
+			// the two-stage 44100->96000 case by 31 (672 versus 703).
+			const tol = 4
+			if got < measured-tol || got > measured+tol {
+				t.Errorf("%v to %v %v: GetLatency()=%d, measured deficit %d",
+					c.in, c.out, preset, got, measured)
+			}
+		}
+	}
+}
diff --git a/resample.go b/resample.go
index c8198cc..428df91 100644
--- a/resample.go
+++ b/resample.go
@@ -31,8 +31,10 @@ type Resampler interface {
 	// MultiFlusher interface instead to drain every channel.
 	Flush() ([]float64, error)
 
-	// GetLatency returns the resampler latency in samples.
-	// This is the delay between input and output due to filtering.
+	// GetLatency returns the startup deficit in output samples: how many
+	// samples early Process calls withhold while internal filters prime.
+	// Real-time users can prime an output FIFO with this many samples of
+	// silence.
 	GetLatency() int
 
 	// Reset clears all internal state and buffers.

From 7e9f7a909452f3b41903ad51111a6d2e444671c6 Mon Sep 17 00:00:00 2001
From: "Tomi P. Hakala" <7030001+tphakala@users.noreply.github.com>
Date: Fri, 17 Jul 2026 20:09:12 +0300
Subject: [PATCH 4/6] refactor: retire duplicate polyphase design path and dead
 LinearStage

internal/filter/polyphase.go carried a parallel polyphase design
implementation that diverged from the engine (edge clamping, DC-gain
normalization, unchecked GetCoefficient, duplicated cubic constants)
and was consumed only by the cmd/analyze-filter diagnostic, so the tool
analyzed a filter the product never ships. Remove both along with the
production-dead LinearStage and its constants; the engine's design code
in internal/engine is now the single polyphase implementation. kaiser.go
stays, the engine depends on it.

Refs #55
---
 .github/workflows/ci.yml                  |   1 -
 CHANGELOG.md                              |   9 +
 cmd/analyze-filter/analyze_filter_gain.go | 132 -------
 codecov.yml                               |   1 -
 internal/engine/buffer_integrity_test.go  |  35 --
 internal/engine/constants.go              |  12 -
 internal/engine/cubic.go                  |  92 -----
 internal/engine/edge_cases_test.go        |   9 -
 internal/engine/extra_engine_test.go      |  17 +-
 internal/engine/reset_state_test.go       |   9 +-
 internal/filter/polyphase.go              | 385 --------------------
 internal/filter/polyphase_test.go         | 413 ----------------------
 12 files changed, 12 insertions(+), 1103 deletions(-)
 delete mode 100644 cmd/analyze-filter/analyze_filter_gain.go
 delete mode 100644 internal/filter/polyphase.go
 delete mode 100644 internal/filter/polyphase_test.go

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index c88713a..eefd001 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -63,7 +63,6 @@ jobs:
         run: |
           go build ./cmd/resample-wav
           go build ./cmd/resample
-          go build ./cmd/analyze-filter
       - name: Verify go.mod is tidy
         run: |
           go mod tidy
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 58391a6..cc79087 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -57,6 +57,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   interpolation (matching `New()` and the documented contract) instead of a
   full FIR pipeline; latency drops accordingly. (#51)
 
+### Removed
+
+- The internal duplicate polyphase filter design path
+  (`internal/filter/polyphase.go`) and the `cmd/analyze-filter` diagnostic
+  that consumed it; the engine's design code in `internal/engine` is the
+  single polyphase implementation. (#55)
+- The production-dead `LinearStage` from `internal/engine`; only its own
+  tests referenced it. (#55)
+
 ## [1.4.0] - 2026-05-29
 
 ### Added
diff --git a/cmd/analyze-filter/analyze_filter_gain.go b/cmd/analyze-filter/analyze_filter_gain.go
deleted file mode 100644
index 24379c5..0000000
--- a/cmd/analyze-filter/analyze_filter_gain.go
+++ /dev/null
@@ -1,132 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Tomi P. Hakala
-// SPDX-License-Identifier: LGPL-2.1-or-later
-
-package main
-
-import (
-	"fmt"
-
-	"github.com/tphakala/go-audio-resampler/internal/filter"
-)
-
-const (
-	// Filter design parameters (matching soxr defaults)
-	defaultNumPhases    = 80    // Number of polyphase filter phases
-	defaultCutoff       = 0.45  // Cutoff frequency for CD→DAT
-	defaultTransitionBW = 0.05  // Transition bandwidth
-	defaultAttenuation  = 100.0 // Stopband attenuation in dB
-
-	// Phase calculation constants
-	phaseShiftBits = 8  // log2(256) = 8 bits for phase indexing
-	fracBits       = 32 // 32-bit fractional precision
-
-	// Display limits
-	maxPhasesToShow = 5    // Maximum phases to display in detail
-	testIterations  = 1000 // Number of iterations for phase usage test
-)
-
-func main() {
-	// Test 2x upsampling case
-	fmt.Println("=== Analyzing Filter DC Gain ===")
-
-	params := filter.PolyphaseParams{
-		NumPhases:    defaultNumPhases,
-		Cutoff:       defaultCutoff,
-		TransitionBW: defaultTransitionBW,
-		Attenuation:  defaultAttenuation,
-		InterpOrder:  filter.InterpLinear,
-		Gain:         1.0,
-	}
-
-	pfb, err := filter.DesignPolyphaseFilterBank(params)
-	if err != nil {
-		fmt.Printf("Error: %v\n", err)
-		return
-	}
-
-	fmt.Printf("Filter bank info:\n")
-	fmt.Printf("  NumPhases: %d\n", pfb.NumPhases)
-	fmt.Printf("  TapsPerPhase: %d\n", pfb.TapsPerPhase)
-	fmt.Printf("  TotalTaps: %d\n", pfb.TotalTaps)
-	fmt.Printf("  InterpOrder: %d\n\n", pfb.InterpOrder)
-
-	// Calculate DC gain of each phase
-	fmt.Println("DC gain per phase:")
-	var totalDC float64
-	phasesToShow := []int{0, 1, 2, 3, 4, 5, 6, 7, 31, 32, 33}
-	phaseGains := make(map[int]float64)
-
-	for phase := range pfb.NumPhases {
-		var phaseDC float64
-		for tap := range pfb.TapsPerPhase {
-			coef := pfb.GetCoefficient(tap, phase, 0.0)
-			phaseDC += coef
-		}
-		phaseGains[phase] = phaseDC
-		totalDC += phaseDC
-	}
-
-	for _, phase := range phasesToShow {
-		if phase < pfb.NumPhases {
-			fmt.Printf("  Phase %2d: %.10f\n", phase, phaseGains[phase])
-		}
-	}
-	if pfb.NumPhases > len(phasesToShow) {
-		fmt.Printf("  ... (%d more phases)\n", pfb.NumPhases-len(phasesToShow))
-	}
-
-	fmt.Printf("\nTotal DC gain (sum of all phases): %.10f\n", totalDC)
-	fmt.Printf("Average DC gain per phase: %.10f\n", totalDC/float64(pfb.NumPhases))
-
-	// Test multiple ratios
-	testRatios := []struct {
-		ratio float64
-		name  string
-	}{
-		{2.0, "2x upsampling"},
-		{0.5, "2x downsampling"},
-		{48000.0 / 44100.0, "CD→DAT"},
-		{44100.0 / 48000.0, "DAT→CD"},
-		{1.5, "3:2 upsampling"},
-	}
-
-	phaseShift := uint(phaseShiftBits)
-	fracScale := float64(uint64(1) << fracBits)
-
-	for _, test := range testRatios {
-		fmt.Printf("\n=== %s (ratio = %.6f) ===\n", test.name, test.ratio)
-
-		// Calculate phase step
-		phaseStep := uint32(fracScale / test.ratio)
-		fmt.Printf("  Phase step: 0x%08x\n", phaseStep)
-
-		// Simulate which phases are used
-		phaseFrac := uint32(0)
-		usedPhases := make(map[int]bool)
-		var sumUsedPhaseDC float64
-
-		// Test enough outputs to find the pattern
-		for range testIterations {
-			phaseIndex := int(phaseFrac >> (fracBits - phaseShift))
-			if phaseIndex >= pfb.NumPhases {
-				phaseIndex = pfb.NumPhases - 1
-			}
-			if !usedPhases[phaseIndex] {
-				usedPhases[phaseIndex] = true
-				sumUsedPhaseDC += phaseGains[phaseIndex]
-				if len(usedPhases) <= maxPhasesToShow {
-					fmt.Printf("    Phase %d: DC gain = %.10f\n",
-						phaseIndex, phaseGains[phaseIndex])
-				}
-			}
-			phaseFrac += phaseStep
-		}
-
-		numUsedPhases := len(usedPhases)
-		avgUsedPhaseDC := sumUsedPhaseDC / float64(numUsedPhases)
-		fmt.Printf("  Used %d unique phases (out of %d)\n", numUsedPhases, pfb.NumPhases)
-		fmt.Printf("  Average DC gain of used phases: %.10f\n", avgUsedPhaseDC)
-		fmt.Printf("  Correct outputGain: 1.0 / %.10f = %.10f\n",
-			avgUsedPhaseDC, 1.0/avgUsedPhaseDC)
-	}
-}
diff --git a/codecov.yml b/codecov.yml
index 64c5096..c7f675e 100644
--- a/codecov.yml
+++ b/codecov.yml
@@ -8,4 +8,3 @@ ignore:
   - "examples/**"
   - "internal/testutil/**"
   - "cmd/resample/**"
-  - "cmd/analyze-filter/**"
diff --git a/internal/engine/buffer_integrity_test.go b/internal/engine/buffer_integrity_test.go
index 03613e1..3d9db64 100644
--- a/internal/engine/buffer_integrity_test.go
+++ b/internal/engine/buffer_integrity_test.go
@@ -438,38 +438,3 @@ func TestCubicStage_BufferIntegrity(t *testing.T) {
 
 	t.Log("Cubic stage buffer integrity verified")
 }
-
-// TestLinearStage_BufferIntegrity verifies LinearStage doesn't have buffer issues.
-func TestLinearStage_BufferIntegrity(t *testing.T) {
-	stage := NewLinearStage(2.0)
-
-	// Generate test signal
-	input := make([]float64, 1000)
-	for i := range input {
-		input[i] = math.Sin(2.0 * math.Pi * float64(i) / 100)
-	}
-
-	// First process call
-	output1, err := stage.Process(input)
-	require.NoError(t, err, "First Process() failed")
-
-	// Save values
-	savedOutput := make([]float64, len(output1))
-	copy(savedOutput, output1)
-
-	// Second process call
-	input2 := make([]float64, 500)
-	for i := range input2 {
-		input2[i] = math.Cos(2.0 * math.Pi * float64(i) / 50)
-	}
-	_, err = stage.Process(input2)
-	require.NoError(t, err, "Second Process() failed")
-
-	// Verify output1 was not corrupted
-	for i, expected := range savedOutput {
-		assert.InDelta(t, expected, output1[i], 1e-15,
-			"output1[%d] was corrupted", i)
-	}
-
-	t.Log("Linear stage buffer integrity verified")
-}
diff --git a/internal/engine/constants.go b/internal/engine/constants.go
index 85d762c..3997cf3 100644
--- a/internal/engine/constants.go
+++ b/internal/engine/constants.go
@@ -14,15 +14,3 @@ const (
 	// Memory usage estimate for cubic stage (bytes)
 	cubicMemoryUsage = 64
 )
-
-// Linear interpolation constants
-const (
-	// Linear interpolation uses 2-point window
-	linearInterpolationPoints = 2
-
-	// Linear interpolation latency
-	linearLatencySamples = 1
-
-	// Memory usage estimate for linear stage (bytes)
-	linearMemoryUsage = 32
-)
diff --git a/internal/engine/cubic.go b/internal/engine/cubic.go
index 0704fc7..86ac396 100644
--- a/internal/engine/cubic.go
+++ b/internal/engine/cubic.go
@@ -212,95 +212,3 @@ func (c *CubicStage[F]) GetPhases() int {
 func (c *CubicStage[F]) GetSIMDInfo() string {
 	return ""
 }
-
-// LinearStage implements linear (2-point, 1st order) interpolation.
-// Even faster than cubic but lower quality.
-type LinearStage struct {
-	ratio   float64
-	phase   float64
-	prev    float64
-	latency int
-}
-
-// NewLinearStage creates a new linear interpolation stage.
-func NewLinearStage(ratio float64) *LinearStage {
-	return &LinearStage{
-		ratio:   ratio,
-		phase:   0,
-		latency: linearLatencySamples,
-	}
-}
-
-// Process resamples input using linear interpolation.
-func (l *LinearStage) Process(input []float64) ([]float64, error) {
-	if len(input) == 0 {
-		return []float64{}, nil
-	}
-
-	outputSize := int(math.Ceil(float64(len(input)) * l.ratio))
-	output := make([]float64, 0, outputSize)
-
-	for _, sample := range input {
-		// Generate output samples between previous and current
-		for l.phase < 1.0 {
-			// Linear interpolation: y = (1-x)*prev + x*current
-			y := (1-l.phase)*l.prev + l.phase*sample
-			output = append(output, y)
-
-			// Advance phase
-			l.phase += 1.0 / l.ratio
-		}
-
-		// Update state
-		l.prev = sample
-		l.phase -= 1.0
-	}
-
-	return output, nil
-}
-
-// Flush returns any remaining samples.
-func (l *LinearStage) Flush() ([]float64, error) {
-	return []float64{}, nil
-}
-
-// Reset clears internal state.
-func (l *LinearStage) Reset() {
-	l.phase = 0
-	l.prev = 0
-}
-
-// GetRatio returns the stage's resampling ratio.
-func (l *LinearStage) GetRatio() float64 {
-	return l.ratio
-}
-
-// GetLatency returns the stage latency in samples.
-func (l *LinearStage) GetLatency() int {
-	return l.latency
-}
-
-// GetMinInput returns the minimum input size for processing.
-func (l *LinearStage) GetMinInput() int {
-	return 1
-}
-
-// GetMemoryUsage returns approximate memory usage in bytes.
-func (l *LinearStage) GetMemoryUsage() int64 {
-	return linearMemoryUsage
-}
-
-// GetFilterLength returns 0 as linear doesn't use a filter.
-func (l *LinearStage) GetFilterLength() int {
-	return linearInterpolationPoints
-}
-
-// GetPhases returns 0 as linear doesn't use phases.
-func (l *LinearStage) GetPhases() int {
-	return 0
-}
-
-// GetSIMDInfo returns empty as linear doesn't use SIMD.
-func (l *LinearStage) GetSIMDInfo() string {
-	return ""
-}
diff --git a/internal/engine/edge_cases_test.go b/internal/engine/edge_cases_test.go
index dc82f8d..0327d70 100644
--- a/internal/engine/edge_cases_test.go
+++ b/internal/engine/edge_cases_test.go
@@ -84,15 +84,6 @@ func TestCubicStage_EmptyInput(t *testing.T) {
 	assert.Empty(t, output, "Output should be empty for empty input")
 }
 
-// TestLinearStage_EmptyInput verifies LinearStage handles empty input correctly.
-func TestLinearStage_EmptyInput(t *testing.T) {
-	stage := NewLinearStage(2.0)
-
-	output, err := stage.Process([]float64{})
-	require.NoError(t, err, "Process() with empty input should not error")
-	assert.Empty(t, output, "Output should be empty for empty input")
-}
-
 // =============================================================================
 // Single Sample Tests
 // =============================================================================
diff --git a/internal/engine/extra_engine_test.go b/internal/engine/extra_engine_test.go
index 057fb95..88497d2 100644
--- a/internal/engine/extra_engine_test.go
+++ b/internal/engine/extra_engine_test.go
@@ -10,7 +10,7 @@ import (
 	"github.com/stretchr/testify/require"
 )
 
-func TestCubicAndLinearStageGetters(t *testing.T) {
+func TestCubicStageGetters(t *testing.T) {
 	// CubicStage float64
 	cs := NewCubicStage[float64](1.5)
 	assert.InDelta(t, 1.5, cs.GetRatio(), 1e-9)
@@ -20,21 +20,6 @@ func TestCubicAndLinearStageGetters(t *testing.T) {
 	assert.Equal(t, 4, cs.GetFilterLength()) // Hermite cubic interpolation uses 4 points
 	assert.Equal(t, 0, cs.GetPhases())
 	assert.Empty(t, cs.GetSIMDInfo())
-
-	// LinearStage
-	ls := NewLinearStage(1.5)
-	assert.InDelta(t, 1.5, ls.GetRatio(), 1e-9)
-	assert.Equal(t, linearLatencySamples, ls.GetLatency())
-	assert.Equal(t, 1, ls.GetMinInput())
-	assert.Equal(t, int64(linearMemoryUsage), ls.GetMemoryUsage())
-	assert.Equal(t, 2, ls.GetFilterLength()) // Linear interpolation uses 2 points
-	assert.Equal(t, 0, ls.GetPhases())
-	assert.Empty(t, ls.GetSIMDInfo())
-
-	// LinearStage Flush
-	flushed, err := ls.Flush()
-	require.NoError(t, err)
-	assert.Empty(t, flushed)
 }
 
 func TestStageAdapterGetters(t *testing.T) {
diff --git a/internal/engine/reset_state_test.go b/internal/engine/reset_state_test.go
index 894ed7d..d89c4d7 100644
--- a/internal/engine/reset_state_test.go
+++ b/internal/engine/reset_state_test.go
@@ -145,13 +145,13 @@ func TestResampler_Reset(t *testing.T) {
 	}
 }
 
-// interpolationStage is an interface for testing cubic and linear stages.
+// interpolationStage is an interface for testing interpolation stages.
 type interpolationStage interface {
 	Process([]float64) ([]float64, error)
 	Reset()
 }
 
-// TestInterpolationStages_Reset verifies Reset() properly clears CubicStage and LinearStage state.
+// TestInterpolationStages_Reset verifies Reset() properly clears CubicStage state.
 func TestInterpolationStages_Reset(t *testing.T) {
 	testCases := []struct {
 		name     string
@@ -163,11 +163,6 @@ func TestInterpolationStages_Reset(t *testing.T) {
 			newStage: func() interpolationStage { return NewCubicStage[float64](2.0) },
 			newFresh: func() interpolationStage { return NewCubicStage[float64](2.0) },
 		},
-		{
-			name:     "LinearStage",
-			newStage: func() interpolationStage { return NewLinearStage(2.0) },
-			newFresh: func() interpolationStage { return NewLinearStage(2.0) },
-		},
 	}
 
 	for _, tc := range testCases {
diff --git a/internal/filter/polyphase.go b/internal/filter/polyphase.go
deleted file mode 100644
index 4281812..0000000
--- a/internal/filter/polyphase.go
+++ /dev/null
@@ -1,385 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Tomi P. Hakala
-// SPDX-License-Identifier: LGPL-2.1-or-later
-
-package filter
-
-import (
-	"fmt"
-	"math"
-
-	"github.com/tphakala/go-audio-resampler/internal/mathutil"
-)
-
-const (
-	// Default number of phases for polyphase filter bank
-	defaultNumPhases = 256
-
-	// Interpolation orders
-	interpOrderNone   = 0 // No coefficient interpolation
-	interpOrderLinear = 1 // Linear interpolation between phases
-	interpOrderCubic  = 3 // Cubic interpolation between phases
-
-	// Coefficient storage multiplier
-	minNumPhases = 2
-	maxNumPhases = 8192
-
-	// Minimum taps per phase for effective filtering
-	// With too few taps, the filter cannot provide good stopband attenuation
-	// and DC gain varies significantly between phases
-	// Increased from 4 to 16 to match soxr's approach and ensure uniform DC gain
-	minTapsPerPhase = 16
-
-	// Polyphase decomposition constants
-	nextPhaseOffset       = 1
-	secondNextPhaseOffset = 2
-
-	// Interpolation polynomial coefficients
-	cubicCenterCoeff = 0.5
-	cubicDCoeff      = 1.0 / 6.0
-	cubicCMultiplier = 4.0
-
-	// Frequency response calculation
-	frequencyNyquistDivisor = 2
-)
-
-// InterpOrder represents the coefficient interpolation order.
-type InterpOrder int
-
-const (
-	// InterpNone means no interpolation (nearest phase)
-	InterpNone InterpOrder = interpOrderNone
-	// InterpLinear means linear interpolation between adjacent phases
-	InterpLinear InterpOrder = interpOrderLinear
-	// InterpCubic means cubic interpolation between phases
-	InterpCubic InterpOrder = interpOrderCubic
-)
-
-// PolyphaseFilterBank represents a polyphase decomposition of an FIR filter.
-//
-// The filter is decomposed into multiple phases, where each phase contains
-// a decimated version of the original filter. This allows efficient arbitrary
-// ratio resampling using fixed-point phase accumulation.
-//
-// Coefficient storage format (per tap, per phase):
-//   - InterpNone:   [coef]
-//   - InterpLinear: [coef, delta]  where value = coef + delta*x
-//   - InterpCubic:  [coef, b, c, d] where value = coef + (b + (c + d*x)*x)*x
-type PolyphaseFilterBank struct {
-	// Coeffs stores all filter coefficients in a flat array.
-	// Layout: [phase0_tap0_coefs...][phase0_tap1_coefs...]...[phaseN_tapM_coefs...]
-	// Each tap stores (InterpOrder+1) coefficients for polynomial evaluation.
-	Coeffs []float64
-
-	// NumPhases is the number of phases (polyphase branches)
-	NumPhases int
-
-	// TapsPerPhase is the number of taps in each phase
-	TapsPerPhase int
-
-	// TotalTaps is the original filter length before decomposition
-	TotalTaps int
-
-	// InterpOrder is the coefficient interpolation order (0, 1, or 3)
-	InterpOrder InterpOrder
-
-	// Cutoff is the normalized cutoff frequency used in design
-	Cutoff float64
-
-	// Attenuation is the stopband attenuation in dB
-	Attenuation float64
-}
-
-// PolyphaseParams holds parameters for polyphase filter bank design.
-type PolyphaseParams struct {
-	// NumPhases is the number of polyphase branches (e.g., 64, 256, 1024)
-	// Higher values allow finer phase resolution and better quality
-	NumPhases int
-
-	// Cutoff is the normalized cutoff frequency (0 to 0.5)
-	// For upsampling: typically 0.5/upsampleRatio
-	// For downsampling: typically 0.5
-	Cutoff float64
-
-	// TransitionBW is the transition bandwidth as fraction of sample rate
-	TransitionBW float64
-
-	// Attenuation is the desired stopband attenuation in dB
-	Attenuation float64
-
-	// InterpOrder specifies coefficient interpolation (0, 1, or 3)
-	// 0: No interpolation (nearest phase)
-	// 1: Linear interpolation (~6dB improvement)
-	// 3: Cubic interpolation (~12dB improvement)
-	InterpOrder InterpOrder
-
-	// Gain is the passband gain (typically 1.0)
-	Gain float64
-}
-
-// Validate checks if polyphase parameters are valid.
-func (pp *PolyphaseParams) Validate() error {
-	if pp.NumPhases < minNumPhases || pp.NumPhases > maxNumPhases {
-		return fmt.Errorf("number of phases %d out of range [%d, %d]",
-			pp.NumPhases, minNumPhases, maxNumPhases)
-	}
-
-	if pp.Cutoff <= 0 || pp.Cutoff >= 0.5 {
-		return fmt.Errorf("cutoff frequency %f out of range (0, 0.5)", pp.Cutoff)
-	}
-
-	if pp.TransitionBW <= 0 || pp.TransitionBW >= 0.5 {
-		return fmt.Errorf("transition bandwidth %f out of range (0, 0.5)", pp.TransitionBW)
-	}
-
-	if pp.Attenuation < 0 {
-		return fmt.Errorf("attenuation %f dB must be positive", pp.Attenuation)
-	}
-
-	if pp.InterpOrder != InterpNone && pp.InterpOrder != InterpLinear && pp.InterpOrder != InterpCubic {
-		return fmt.Errorf("invalid interpolation order %d (must be 0, 1, or 3)", pp.InterpOrder)
-	}
-
-	if pp.Gain <= 0 {
-		return fmt.Errorf("gain %f must be positive", pp.Gain)
-	}
-
-	return nil
-}
-
-// DesignPolyphaseFilterBank creates a polyphase filter bank from the given parameters.
-//
-// The process:
-// 1. Design a prototype lowpass filter using Kaiser window method
-// 2. Decompose the filter into multiple phases
-// 3. Compute interpolation coefficients for sub-phase precision
-//
-// Returns the polyphase filter bank ready for use in resampling.
-func DesignPolyphaseFilterBank(params PolyphaseParams) (*PolyphaseFilterBank, error) {
-	if err := params.Validate(); err != nil {
-		return nil, fmt.Errorf("invalid polyphase parameters: %w", err)
-	}
-
-	// Calculate minimum filter length to ensure adequate taps per phase
-	minTotalTaps := minTapsPerPhase * params.NumPhases
-
-	// Design prototype lowpass filter with adequate length
-	// Use the auto design first to get estimated length
-	estimatedFilter, err := DesignLowPassFilterAuto(
-		params.Cutoff,
-		params.TransitionBW,
-		params.Attenuation,
-		params.Gain,
-	)
-	if err != nil {
-		return nil, fmt.Errorf("failed to design prototype filter: %w", err)
-	}
-
-	// If the estimated filter is too short, design a longer one by using the minimum length
-	var prototypeFilter []float64
-	if len(estimatedFilter) < minTotalTaps {
-		// Design with explicit length
-		beta := mathutil.KaiserBeta(params.Attenuation)
-		window := KaiserWindow(minTotalTaps, beta)
-
-		// Generate windowed sinc
-		prototypeFilter = make([]float64, minTotalTaps)
-		center := float64(minTotalTaps-1) / windowNormalizationFactor
-
-		for n := range minTotalTaps {
-			x := float64(n) - center
-			var sincValue float64
-			const sincZeroThreshold = 1e-10
-			if math.Abs(x) < sincZeroThreshold {
-				sincValue = windowNormalizationFactor * params.Cutoff
-			} else {
-				arg := windowNormalizationFactor * math.Pi * params.Cutoff * x
-				sincValue = math.Sin(arg) / (math.Pi * x)
-			}
-			prototypeFilter[n] = sincValue * window[n]
-		}
-
-		// Normalize for desired gain
-		// Scale by NumPhases so that the average DC gain per phase is 1.0
-		// This matches soxr's approach: the prototype filter has total DC gain = NumPhases
-		sum := 0.0
-		for _, coeff := range prototypeFilter {
-			sum += coeff
-		}
-		if math.Abs(sum) > sincZeroThreshold {
-			scale := params.Gain * float64(params.NumPhases) / sum
-			for i := range prototypeFilter {
-				prototypeFilter[i] *= scale
-			}
-		}
-	} else {
-		prototypeFilter = estimatedFilter
-	}
-
-	// Create polyphase filter bank
-	pfb := &PolyphaseFilterBank{
-		NumPhases:   params.NumPhases,
-		TotalTaps:   len(prototypeFilter),
-		InterpOrder: params.InterpOrder,
-		Cutoff:      params.Cutoff,
-		Attenuation: params.Attenuation,
-	}
-
-	// Calculate taps per phase
-	pfb.TapsPerPhase = (pfb.TotalTaps + pfb.NumPhases - 1) / pfb.NumPhases
-
-	// Decompose prototype filter into phases with interpolation coefficients
-	pfb.Coeffs = decomposePolyphase(prototypeFilter, pfb.NumPhases, pfb.TapsPerPhase, pfb.InterpOrder)
-
-	return pfb, nil
-}
-
-// decomposePolyphase decomposes a prototype filter into polyphase branches
-// and computes interpolation coefficients.
-//
-// The polyphase decomposition distributes the prototype filter coefficients
-// across multiple phases. For coefficient interpolation, we compute polynomial
-// coefficients that allow smooth interpolation between adjacent phases.
-func decomposePolyphase(prototype []float64, numPhases, tapsPerPhase int, interpOrder InterpOrder) []float64 {
-	// Allocate coefficient storage
-	// Each tap in each phase stores (interpOrder+1) coefficients
-	coeffsPerTap := int(interpOrder) + 1
-	totalCoeffs := tapsPerPhase * numPhases * coeffsPerTap
-	coeffs := make([]float64, totalCoeffs)
-
-	// Helper function to get prototype coefficient with boundary handling
-	getProtoCoeff := func(phase, tap int) float64 {
-		idx := tap*numPhases + phase
-		if idx < 0 || idx >= len(prototype) {
-			return 0.0
-		}
-		return prototype[idx]
-	}
-
-	// For each tap position in each phase, compute interpolation coefficients
-	for tap := range tapsPerPhase {
-		for phase := range numPhases {
-			// Get coefficients from adjacent phases for interpolation
-			// f0 = current phase, f1 = next phase, etc.
-			prevPhase := max(phase-1, 0)
-
-			f0 := getProtoCoeff(phase, tap)
-			f1 := getProtoCoeff(phase+nextPhaseOffset, tap)
-			fm1 := getProtoCoeff(prevPhase, tap)
-			f2 := getProtoCoeff(phase+secondNextPhaseOffset, tap)
-
-			// Calculate base index for this tap/phase combination
-			baseIdx := (tap*numPhases + phase) * coeffsPerTap
-
-			// Store coefficients based on interpolation order
-			switch interpOrder {
-			case InterpNone:
-				// No interpolation: just store the coefficient
-				coeffs[baseIdx] = f0
-
-			case InterpLinear:
-				// Linear interpolation: f(x) = f0 + b*x
-				// where b = f1 - f0
-				coeffs[baseIdx] = f0
-				coeffs[baseIdx+1] = f1 - f0
-
-			case InterpCubic:
-				// Cubic interpolation: f(x) = f0 + b*x + c*x^2 + d*x^3
-				// Using centered finite differences for smooth interpolation
-				c := cubicCenterCoeff*(f1+fm1) - f0
-				d := cubicDCoeff * (f2 - f1 + fm1 - f0 - cubicCMultiplier*c)
-				b := f1 - f0 - d - c
-
-				coeffs[baseIdx] = f0
-				coeffs[baseIdx+1] = b
-				coeffs[baseIdx+2] = c
-				coeffs[baseIdx+3] = d
-			}
-		}
-	}
-
-	return coeffs
-}
-
-// GetCoefficient returns the interpolated coefficient for a given tap and fractional phase.
-//
-// Parameters:
-//   - tap: The tap index (0 to TapsPerPhase-1)
-//   - phase: The integer phase index (0 to NumPhases-1)
-//   - frac: The fractional phase position [0, 1) for sub-phase interpolation
-func (pfb *PolyphaseFilterBank) GetCoefficient(tap, phase int, frac float64) float64 {
-	coeffsPerTap := int(pfb.InterpOrder) + 1
-	baseIdx := (tap*pfb.NumPhases + phase) * coeffsPerTap
-
-	switch pfb.InterpOrder {
-	case InterpNone:
-		return pfb.Coeffs[baseIdx]
-
-	case InterpLinear:
-		// Linear: f0 + b*x
-		f0 := pfb.Coeffs[baseIdx]
-		b := pfb.Coeffs[baseIdx+1]
-		return f0 + b*frac
-
-	case InterpCubic:
-		// Cubic: f0 + (b + (c + d*x)*x)*x
-		// Horner's method for efficient evaluation
-		f0 := pfb.Coeffs[baseIdx]
-		b := pfb.Coeffs[baseIdx+1]
-		c := pfb.Coeffs[baseIdx+2]
-		d := pfb.Coeffs[baseIdx+3]
-		return f0 + (b+(c+d*frac)*frac)*frac
-
-	default:
-		return pfb.Coeffs[baseIdx]
-	}
-}
-
-// ComputeFrequencyResponse computes the frequency response of the polyphase filter bank.
-// This evaluates the response of a single phase (phase 0) as a representative.
-func (pfb *PolyphaseFilterBank) ComputeFrequencyResponse(numPoints int) FilterResponse {
-	if numPoints <= 0 {
-		numPoints = 512
-	}
-
-	response := FilterResponse{
-		Frequencies: make([]float64, numPoints),
-		Magnitude:   make([]float64, numPoints),
-		Phase:       make([]float64, numPoints),
-	}
-
-	// Extract phase 0 coefficients for frequency response calculation
-	phase0Coeffs := make([]float64, pfb.TapsPerPhase)
-	coeffsPerTap := int(pfb.InterpOrder) + 1
-
-	for tap := range pfb.TapsPerPhase {
-		baseIdx := (tap*pfb.NumPhases + 0) * coeffsPerTap
-		phase0Coeffs[tap] = pfb.Coeffs[baseIdx]
-	}
-
-	// Compute DTFT of phase 0
-	const twoPi = 2.0 * math.Pi
-	for k := range numPoints {
-		freq := float64(k) / float64(frequencyNyquistDivisor*numPoints)
-		response.Frequencies[k] = freq
-
-		var realPart, imagPart float64
-		omega := twoPi * freq
-
-		for n, h := range phase0Coeffs {
-			angle := omega * float64(n*pfb.NumPhases)
-			realPart += h * math.Cos(angle)
-			imagPart -= h * math.Sin(angle)
-		}
-
-		response.Magnitude[k] = math.Sqrt(realPart*realPart + imagPart*imagPart)
-		response.Phase[k] = math.Atan2(imagPart, realPart)
-	}
-
-	return response
-}
-
-// GetMemoryUsage returns the approximate memory usage in bytes.
-func (pfb *PolyphaseFilterBank) GetMemoryUsage() int64 {
-	const bytesPerFloat64 = 8
-	return int64(len(pfb.Coeffs)) * bytesPerFloat64
-}
diff --git a/internal/filter/polyphase_test.go b/internal/filter/polyphase_test.go
deleted file mode 100644
index 8932a54..0000000
--- a/internal/filter/polyphase_test.go
+++ /dev/null
@@ -1,413 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Tomi P. Hakala
-// SPDX-License-Identifier: LGPL-2.1-or-later
-
-package filter
-
-import (
-	"fmt"
-	"math"
-	"testing"
-
-	"github.com/stretchr/testify/assert"
-	"github.com/stretchr/testify/require"
-	"github.com/tphakala/go-audio-resampler/internal/testutil"
-)
-
-const (
-	// Test parameters for polyphase tests
-	testNumPhases64   = 64
-	testNumPhases256  = 256
-	testNumPhases1024 = 1024
-
-	testTransition005 = 0.05
-)
-
-// TestPolyphaseParams_Validate tests parameter validation.
-func TestPolyphaseParams_Validate(t *testing.T) {
-	tests := []struct {
-		name    string
-		params  PolyphaseParams
-		wantErr bool
-	}{
-		{
-			name: "valid_params",
-			params: PolyphaseParams{
-				NumPhases:    testNumPhases256,
-				Cutoff:       testCutoff0_25,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  InterpLinear,
-				Gain:         testGainUnity,
-			},
-			wantErr: false,
-		},
-		{
-			name: "too_few_phases",
-			params: PolyphaseParams{
-				NumPhases:    1,
-				Cutoff:       testCutoff0_25,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  InterpLinear,
-				Gain:         testGainUnity,
-			},
-			wantErr: true,
-		},
-		{
-			name: "too_many_phases",
-			params: PolyphaseParams{
-				NumPhases:    10000,
-				Cutoff:       testCutoff0_25,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  InterpLinear,
-				Gain:         testGainUnity,
-			},
-			wantErr: true,
-		},
-		{
-			name: "invalid_cutoff_low",
-			params: PolyphaseParams{
-				NumPhases:    testNumPhases256,
-				Cutoff:       0.0,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  InterpLinear,
-				Gain:         testGainUnity,
-			},
-			wantErr: true,
-		},
-		{
-			name: "invalid_cutoff_high",
-			params: PolyphaseParams{
-				NumPhases:    testNumPhases256,
-				Cutoff:       0.5,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  InterpLinear,
-				Gain:         testGainUnity,
-			},
-			wantErr: true,
-		},
-		{
-			name: "invalid_interp_order",
-			params: PolyphaseParams{
-				NumPhases:    testNumPhases256,
-				Cutoff:       testCutoff0_25,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  InterpOrder(2), // Invalid: only 0, 1, 3 allowed
-				Gain:         testGainUnity,
-			},
-			wantErr: true,
-		},
-	}
-
-	for _, tt := range tests {
-		t.Run(tt.name, func(t *testing.T) {
-			err := tt.params.Validate()
-			if tt.wantErr {
-				assert.Error(t, err, "expected validation error")
-			} else {
-				assert.NoError(t, err, "unexpected validation error")
-			}
-		})
-	}
-}
-
-// TestDesignPolyphaseFilterBank tests basic polyphase filter bank design.
-func TestDesignPolyphaseFilterBank(t *testing.T) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation80,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	pfb, err := DesignPolyphaseFilterBank(params)
-	require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-	// Check basic properties
-	assert.Equal(t, params.NumPhases, pfb.NumPhases, "NumPhases mismatch")
-	assert.Equal(t, params.InterpOrder, pfb.InterpOrder, "InterpOrder mismatch")
-	assert.Positive(t, pfb.TotalTaps, "TotalTaps should be > 0")
-	assert.Positive(t, pfb.TapsPerPhase, "TapsPerPhase should be > 0")
-
-	// Check coefficient storage size
-	expectedCoeffs := pfb.TapsPerPhase * pfb.NumPhases * (int(pfb.InterpOrder) + 1)
-	assert.Len(t, pfb.Coeffs, expectedCoeffs, "Coeffs length mismatch")
-}
-
-// TestPolyphaseFilterBank_InterpolationOrders tests all interpolation orders.
-func TestPolyphaseFilterBank_InterpolationOrders(t *testing.T) {
-	orders := []struct {
-		name  string
-		order InterpOrder
-	}{
-		{"none", InterpNone},
-		{"linear", InterpLinear},
-		{"cubic", InterpCubic},
-	}
-
-	for _, ord := range orders {
-		t.Run(ord.name, func(t *testing.T) {
-			params := PolyphaseParams{
-				NumPhases:    testNumPhases64,
-				Cutoff:       testCutoff0_25,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  ord.order,
-				Gain:         testGainUnity,
-			}
-
-			pfb, err := DesignPolyphaseFilterBank(params)
-			require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-			// Verify coefficient storage
-			coeffsPerTap := int(ord.order) + 1
-			expectedSize := pfb.TapsPerPhase * pfb.NumPhases * coeffsPerTap
-			assert.Len(t, pfb.Coeffs, expectedSize, "Coeffs length mismatch")
-		})
-	}
-}
-
-// TestPolyphaseFilterBank_GetCoefficient tests coefficient retrieval and interpolation.
-func TestPolyphaseFilterBank_GetCoefficient(t *testing.T) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation80,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	pfb, err := DesignPolyphaseFilterBank(params)
-	require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-	// Test coefficient retrieval
-	tap := 0
-	phase := 0
-
-	// At frac = 0, should return base coefficient
-	coef0 := pfb.GetCoefficient(tap, phase, 0.0)
-	assert.False(t, math.IsNaN(coef0), "GetCoefficient returned NaN for frac=0")
-
-	// At frac = 1, should approach next phase
-	coef1 := pfb.GetCoefficient(tap, phase, 1.0)
-	assert.False(t, math.IsNaN(coef1), "GetCoefficient returned NaN for frac=1")
-
-	// With linear interpolation, frac=0.5 should be average of endpoints
-	coef0_5 := pfb.GetCoefficient(tap, phase, 0.5)
-	assert.False(t, math.IsNaN(coef0_5), "GetCoefficient returned NaN for frac=0.5")
-
-	// For linear interpolation, midpoint should satisfy interpolation property
-	// This is a basic sanity check, not a strict equality test
-	maxEndpoint := math.Max(math.Abs(coef0), math.Abs(coef1)) * 2
-	assert.LessOrEqual(t, math.Abs(coef0_5), maxEndpoint,
-		"Interpolated coefficient suspiciously large: %f (endpoints: %f, %f)",
-		coef0_5, coef0, coef1)
-}
-
-// TestPolyphaseFilterBank_Structure tests that the filter bank has valid structure.
-// Note: In polyphase decomposition with soxr-style scaling, individual tap coefficients
-// across phases may vary significantly (that's normal). We test structural validity instead.
-func TestPolyphaseFilterBank_Structure(t *testing.T) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases64,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation80,
-		InterpOrder:  InterpNone,
-		Gain:         testGainUnity,
-	}
-
-	pfb, err := DesignPolyphaseFilterBank(params)
-	require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-	// Verify structure is consistent
-	assert.Equal(t, params.NumPhases, pfb.NumPhases, "NumPhases mismatch")
-
-	// Verify minimum taps per phase
-	assert.GreaterOrEqual(t, pfb.TapsPerPhase, 2, "TapsPerPhase should be at least 2")
-
-	// Verify coefficient storage size
-	expectedCoeffs := pfb.TapsPerPhase * pfb.NumPhases * (int(pfb.InterpOrder) + 1)
-	assert.Len(t, pfb.Coeffs, expectedCoeffs, "Coeffs length mismatch")
-
-	// Verify all coefficients are valid (not NaN or Inf)
-	testutil.AssertNoNaNOrInf(t, pfb.Coeffs)
-
-	t.Logf("Filter bank: %d phases, %d taps/phase, %d total coefficients",
-		pfb.NumPhases, pfb.TapsPerPhase, len(pfb.Coeffs))
-}
-
-// TestPolyphaseFilterBank_DCGain tests that DC gain is preserved per phase.
-// soxr-style: each phase should have DC gain ≈ 1.0 (not 1/numPhases).
-// This ensures DC preservation regardless of which phases are used during resampling.
-func TestPolyphaseFilterBank_DCGain(t *testing.T) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation80,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	pfb, err := DesignPolyphaseFilterBank(params)
-	require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-	// Check DC gain across multiple phases (soxr-style: each phase ≈ 1.0)
-	// The average DC gain across all phases should be approximately 1.0
-	var totalDCGain float64
-	for phase := range pfb.NumPhases {
-		var phaseSum float64
-		for tap := range pfb.TapsPerPhase {
-			phaseSum += pfb.GetCoefficient(tap, phase, 0.0)
-		}
-		totalDCGain += phaseSum
-	}
-
-	avgDCGain := totalDCGain / float64(pfb.NumPhases)
-	expectedGain := params.Gain
-
-	// Allow some tolerance for filter design imprecision
-	tolerance := 0.5 // Wider tolerance for average across all phases
-	assert.InDelta(t, expectedGain, avgDCGain, tolerance,
-		"Average DC gain mismatch")
-	t.Logf("Average DC gain across %d phases: %.6f", pfb.NumPhases, avgDCGain)
-}
-
-// TestPolyphaseFilterBank_DifferentPhases tests various phase counts.
-func TestPolyphaseFilterBank_DifferentPhases(t *testing.T) {
-	phaseCounts := []int{testNumPhases64, testNumPhases256, testNumPhases1024}
-
-	for _, numPhases := range phaseCounts {
-		t.Run(fmt.Sprintf("phases_%d", numPhases), func(t *testing.T) {
-			params := PolyphaseParams{
-				NumPhases:    numPhases,
-				Cutoff:       testCutoff0_25,
-				TransitionBW: testTransition005,
-				Attenuation:  testAttenuation80,
-				InterpOrder:  InterpLinear,
-				Gain:         testGainUnity,
-			}
-
-			pfb, err := DesignPolyphaseFilterBank(params)
-			require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-			assert.Equal(t, numPhases, pfb.NumPhases, "NumPhases mismatch")
-
-			// Verify we have coefficients for all phases
-			coeffsPerTap := int(pfb.InterpOrder) + 1
-			expectedSize := pfb.TapsPerPhase * numPhases * coeffsPerTap
-			assert.Len(t, pfb.Coeffs, expectedSize, "Coeffs length mismatch")
-		})
-	}
-}
-
-// TestPolyphaseFilterBank_FrequencyResponse tests frequency response computation.
-func TestPolyphaseFilterBank_FrequencyResponse(t *testing.T) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation80,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	pfb, err := DesignPolyphaseFilterBank(params)
-	require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-	response := pfb.ComputeFrequencyResponse(testNumPoints512)
-
-	assert.Len(t, response.Frequencies, testNumPoints512, "response length mismatch")
-
-	// Check that DC response is reasonable
-	dcMagnitude := response.Magnitude[0]
-	assert.Positive(t, dcMagnitude, "DC magnitude should be > 0")
-	assert.False(t, math.IsNaN(dcMagnitude), "DC magnitude should not be NaN")
-
-	// Check that frequencies are in expected range [0, 0.5]
-	testutil.AssertAllInRange(t, response.Frequencies, 0, 0.5)
-}
-
-// TestPolyphaseFilterBank_MemoryUsage tests memory usage calculation.
-func TestPolyphaseFilterBank_MemoryUsage(t *testing.T) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation80,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	pfb, err := DesignPolyphaseFilterBank(params)
-	require.NoError(t, err, "DesignPolyphaseFilterBank failed")
-
-	memUsage := pfb.GetMemoryUsage()
-	const bytesPerFloat64 = 8
-	expectedUsage := int64(len(pfb.Coeffs)) * bytesPerFloat64
-
-	assert.Equal(t, expectedUsage, memUsage, "GetMemoryUsage mismatch")
-}
-
-// BenchmarkDesignPolyphaseFilterBank benchmarks filter bank design.
-func BenchmarkDesignPolyphaseFilterBank(b *testing.B) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation100,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	b.ResetTimer()
-	for b.Loop() {
-		_, _ = DesignPolyphaseFilterBank(params)
-	}
-}
-
-// BenchmarkPolyphaseGetCoefficient benchmarks coefficient retrieval.
-func BenchmarkPolyphaseGetCoefficient(b *testing.B) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation100,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	pfb, _ := DesignPolyphaseFilterBank(params)
-
-	b.ResetTimer()
-	for b.Loop() {
-		_ = pfb.GetCoefficient(0, 0, 0.5)
-	}
-}
-
-// BenchmarkPolyphaseFrequencyResponse benchmarks frequency response computation.
-func BenchmarkPolyphaseFrequencyResponse(b *testing.B) {
-	params := PolyphaseParams{
-		NumPhases:    testNumPhases256,
-		Cutoff:       testCutoff0_25,
-		TransitionBW: testTransition005,
-		Attenuation:  testAttenuation100,
-		InterpOrder:  InterpLinear,
-		Gain:         testGainUnity,
-	}
-
-	pfb, _ := DesignPolyphaseFilterBank(params)
-
-	b.ResetTimer()
-	for b.Loop() {
-		_ = pfb.ComputeFrequencyResponse(testNumPoints512)
-	}
-}

From e84369d16813188be7b6349d392b3e5171b5b04e Mon Sep 17 00:00:00 2001
From: "Tomi P. Hakala" <7030001+tphakala@users.noreply.github.com>
Date: Fri, 17 Jul 2026 20:44:46 +0300
Subject: [PATCH 5/6] refactor: harden startup-deficit contract and align docs
 from review

Gate review follow-ups: README no longer describes the retired
group-delay semantics of GetLatency on the New(config) path; the
duck-typed deficit assertion uses a named startupDeficitStage interface
with compile-time assertions in stages.go so a refactor of the
StageAdapter embedding fails the build instead of silently degrading to
the group-delay fallback; the cubic deficit formula lives only in
CubicStage.StartupDeficit with the engine delegating to it; the config
path latency test now covers QualityQuick and a new test pins the
fallback branch including its downstream-ratio conversion;
TestCubicStageGetters asserts the new getter; stale or imprecise
comments corrected (StageAdapter.GetLatency consumer note, THD margin
arithmetic, duplicated bank-build description, what the bank-equality
pin does and does not guard).

Refs #52 #53 #54
---
 README.md                                  |  2 +-
 constant.go                                | 10 ++++++++-
 internal/engine/cubic.go                   |  4 +++-
 internal/engine/extra_engine_test.go       |  1 +
 internal/engine/phase_wrap_measure_test.go |  6 +++++-
 internal/engine/polyphase_stage.go         |  5 ++---
 internal/engine/quality_regression_test.go |  8 ++++----
 internal/engine/resampler.go               |  2 +-
 internal/engine/stage_adapter.go           |  9 ++++----
 latency_test.go                            | 24 +++++++++++++++++++++-
 stages.go                                  |  8 ++++++++
 11 files changed, 62 insertions(+), 17 deletions(-)

diff --git a/README.md b/README.md
index 210c7ad..7c531e7 100644
--- a/README.md
+++ b/README.md
@@ -148,7 +148,7 @@ func main() {
 
 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.
+`Latency()` is available on `SimpleResampler` and `SimpleResamplerFloat32` (the `NewEngine`/`NewEngineFloat32` path). Resamplers built from `New(config)` expose the same figure through `GetLatency()` (and `GetInfo().Latency`): the startup deficit in output samples, accurate to within a few samples across the pipeline's stages, so it primes a FIFO the same way.
 
 ```go
 r, err := resampling.NewEngineFloat32(44100, 48000, resampling.QualityHigh)
diff --git a/constant.go b/constant.go
index e8865fa..ba2d910 100644
--- a/constant.go
+++ b/constant.go
@@ -404,6 +404,14 @@ func (r *constantRateResampler) FlushMulti() ([][]float64, error) {
 	return output, nil
 }
 
+// startupDeficitStage is the accounting contract stages provide for accurate
+// latency reporting: the un-rounded startup deficit in the stage's own
+// output-sample domain. Compile-time assertions in stages.go keep every
+// production stage type on this path.
+type startupDeficitStage interface {
+	StartupDeficit() float64
+}
+
 // GetLatency returns the pipeline's startup deficit in output samples: how
 // many samples early Process calls withhold while the stage filters prime.
 // Each stage's deficit is converted into the final output rate domain
@@ -419,7 +427,7 @@ func (r *constantRateResampler) GetLatency() int {
 	total := 0.0
 	for i, stage := range ch.stages {
 		var deficit float64
-		if s, ok := stage.(interface{ StartupDeficit() float64 }); ok {
+		if s, ok := stage.(startupDeficitStage); ok {
 			deficit = s.StartupDeficit()
 		} else {
 			// Fallback for stages without deficit accounting: group-delay
diff --git a/internal/engine/cubic.go b/internal/engine/cubic.go
index 86ac396..45b95ed 100644
--- a/internal/engine/cubic.go
+++ b/internal/engine/cubic.go
@@ -193,7 +193,9 @@ func (c *CubicStage[F]) GetMemoryUsage() int64 {
 }
 
 // StartupDeficit returns the startup deficit in output samples, un-rounded.
-// Mirrors the cubic branch of Resampler.Latency for direct pipeline use.
+// This is the single source of the cubic deficit formula: the engine
+// Resampler.StartupDeficit delegates here when its cubic stage is active,
+// and the New(config) pipeline consumes it directly per stage.
 func (c *CubicStage[F]) StartupDeficit() float64 {
 	return cubicLatencySamples * c.ratio
 }
diff --git a/internal/engine/extra_engine_test.go b/internal/engine/extra_engine_test.go
index 88497d2..e5aae0e 100644
--- a/internal/engine/extra_engine_test.go
+++ b/internal/engine/extra_engine_test.go
@@ -15,6 +15,7 @@ func TestCubicStageGetters(t *testing.T) {
 	cs := NewCubicStage[float64](1.5)
 	assert.InDelta(t, 1.5, cs.GetRatio(), 1e-9)
 	assert.Equal(t, cubicLatencySamples, cs.GetLatency())
+	assert.InDelta(t, cubicLatencySamples*1.5, cs.StartupDeficit(), 1e-9)
 	assert.Equal(t, 1, cs.GetMinInput())
 	assert.Equal(t, int64(cubicMemoryUsage), cs.GetMemoryUsage())
 	assert.Equal(t, 4, cs.GetFilterLength()) // Hermite cubic interpolation uses 4 points
diff --git a/internal/engine/phase_wrap_measure_test.go b/internal/engine/phase_wrap_measure_test.go
index df30903..5fd0e47 100644
--- a/internal/engine/phase_wrap_measure_test.go
+++ b/internal/engine/phase_wrap_measure_test.go
@@ -86,7 +86,11 @@ type phaseWrapResult struct {
 // flat=false it reproduces the original wrap: the neighbour phase is kept
 // within the same tap (phase % numPhases). Everything else, the Catmull-Rom
 // cubic math and the reversed tap storage, is identical (shared code), so any
-// output difference is attributable solely to the boundary indexing.
+// output difference is attributable solely to the boundary indexing. Because
+// both sides share buildCubicInterpBanks, the bank-equality assertion pins
+// the boundary-indexing policy only; the Catmull-Rom math itself is guarded
+// by the THD bounds below and by the active-interpolation THD pins in
+// quality_regression_test.go.
 func buildPhaseInterpBanks(fb *polyphaseFilter, numPhases int, flat bool) (a, b, c, d [][]float64) {
 	coeffs := fb.coeffs
 
diff --git a/internal/engine/polyphase_stage.go b/internal/engine/polyphase_stage.go
index 5c8d33b..651ca17 100644
--- a/internal/engine/polyphase_stage.go
+++ b/internal/engine/polyphase_stage.go
@@ -173,9 +173,8 @@ func NewPolyphaseStage[F simdops.Float](ratio, totalIORatio float64, hasPreStage
 		return filterBank.coeffs[idx]
 	}
 
-	// Allocate coefficient arrays with cubic interpolation support
-	// polyCoeffs = a (base), polyCoeffsB = b (linear), polyCoeffsC = c (quadratic), polyCoeffsD = d (cubic)
-	// Interpolation formula: coef(x) = a + x*(b + x*(c + x*d)) where x ∈ [0, 1)
+	// polyCoeffs = a (base), polyCoeffsB = b (linear), polyCoeffsC = c (quadratic), polyCoeffsD = d (cubic);
+	// evaluated as coef(x) = a + x*(b + x*(c + x*d)) where x ∈ [0, 1)
 	polyCoeffs, polyCoeffsB, polyCoeffsC, polyCoeffsD := buildCubicInterpBanks[F](numPhases, tapsPerPhase, getCoeff)
 
 	return &PolyphaseStage[F]{
diff --git a/internal/engine/quality_regression_test.go b/internal/engine/quality_regression_test.go
index 54d36bc..ff0bb11 100644
--- a/internal/engine/quality_regression_test.go
+++ b/internal/engine/quality_regression_test.go
@@ -47,10 +47,10 @@ const (
 	// Active-interpolation ratio THD limits (#54). Measured on 2026-07-17:
 	// 32000->44100 High: -162.05 dB, Medium: -134.95 dB, Low: -144.22 dB;
 	// 44100->64000 High: -151.63 dB, Medium: -135.03 dB, Low: -139.47 dB.
-	// Limits are measured values rounded up with ~5 dB margin. Medium and Low
-	// reuse the regressionMaxTHD_Medium and regressionMaxTHD_Low constants
-	// above since the measured values here fall within margin of those
-	// existing constants; only High needed a dedicated, tighter limit.
+	// High gets a dedicated limit 6.6 dB above its worst measurement. Medium
+	// and Low reuse the regressionMaxTHD_Medium and regressionMaxTHD_Low
+	// constants above, which sit 6.0 and 9.5 dB above the worst measurements
+	// here, so no dedicated constants are needed for them.
 	regressionMaxTHDInterp_High = -145.0 // dB (worst measured: -151.63 dB)
 
 	// Minimum SNR at 1kHz test frequency (more positive = better)
diff --git a/internal/engine/resampler.go b/internal/engine/resampler.go
index 5c36d94..6641046 100644
--- a/internal/engine/resampler.go
+++ b/internal/engine/resampler.go
@@ -372,7 +372,7 @@ func (r *Resampler[F]) GetStatistics() map[string]int64 {
 // wraps this with math.Ceil.
 func (r *Resampler[F]) StartupDeficit() float64 {
 	if r.cubicStage != nil {
-		return cubicLatencySamples * r.ratio
+		return r.cubicStage.StartupDeficit()
 	}
 	deficitIn := 0.0
 	if r.preStage != nil && r.preStage.factor > 1 {
diff --git a/internal/engine/stage_adapter.go b/internal/engine/stage_adapter.go
index 0ba9580..ac3aab3 100644
--- a/internal/engine/stage_adapter.go
+++ b/internal/engine/stage_adapter.go
@@ -41,10 +41,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.
+// This is a filter group-delay heuristic in the input domain, kept for
+// stage-level diagnostics. Latency reporting flows through StartupDeficit
+// (promoted from the embedded Resampler); streaming users priming a FIFO
+// with silence should use the engine Resampler.Latency() accessor, which
+// reports the startup deficit in output samples.
 func (s *StageAdapter[F]) GetLatency() int {
 	latency := 0
 
diff --git a/latency_test.go b/latency_test.go
index 9782c0d..0e99291 100644
--- a/latency_test.go
+++ b/latency_test.go
@@ -78,7 +78,7 @@ func TestGetLatency_ConfigPath_MatchesMeasuredDeficit(t *testing.T) {
 		{44100, 96000},
 		{48000, 44100},
 	} {
-		for _, preset := range []QualityPreset{QualityLow, QualityMedium, QualityHigh} {
+		for _, preset := range []QualityPreset{QualityQuick, QualityLow, QualityMedium, QualityHigh} {
 			cfg := Config{
 				InputRate:  c.in,
 				OutputRate: c.out,
@@ -113,3 +113,25 @@ func TestGetLatency_ConfigPath_MatchesMeasuredDeficit(t *testing.T) {
 		}
 	}
 }
+
+// The duck-typed StartupDeficit fast path covers every production stage type
+// (enforced at compile time in stages.go); the fallback in GetLatency keeps
+// working for stage types without deficit accounting. Pinned here with the
+// test-only stubStage, including the downstream-ratio conversion.
+func TestGetLatency_FallbackWithoutStartupDeficit(t *testing.T) {
+	r := &constantRateResampler{
+		pipeline: &Pipeline{},
+		channels: []*channelResampler{{
+			stages: []Stage{
+				&stubStage{ratio: 2.0, filterLength: 32},
+				&stubStage{ratio: 0.5, filterLength: 8},
+			},
+		}},
+	}
+	// Per-stage fallback deficits in each stage's own output domain:
+	// 32/2*2.0 = 32 and 8/2*0.5 = 2. The first converts through the
+	// downstream ratio 0.5 to 16, so the total is 18.
+	if got := r.GetLatency(); got != 18 {
+		t.Errorf("GetLatency()=%d, want 18", got)
+	}
+}
diff --git a/stages.go b/stages.go
index 69cfd77..b7ba0d7 100644
--- a/stages.go
+++ b/stages.go
@@ -193,4 +193,12 @@ var (
 	_ pipeline.Stage = (*engine.CubicStage[float64])(nil)
 	_ pipeline.Stage = (*engine.StageAdapter[float64])(nil)
 	_ pipeline.Stage = (*stubStage)(nil)
+
+	// Every production stage must provide startup-deficit accounting so
+	// GetLatency never degrades to its group-delay fallback. StageAdapter
+	// satisfies this through its embedded *engine.Resampler; without these
+	// assertions a refactor of that embedding would break latency reporting
+	// silently instead of failing the build.
+	_ startupDeficitStage = (*engine.CubicStage[float64])(nil)
+	_ startupDeficitStage = (*engine.StageAdapter[float64])(nil)
 )

From 63852b42b95afc9f2f4fef1444b708369c03b8b7 Mon Sep 17 00:00:00 2001
From: "Tomi P. Hakala" <7030001+tphakala@users.noreply.github.com>
Date: Fri, 17 Jul 2026 20:53:22 +0300
Subject: [PATCH 6/6] docs(test): attribute phase-boundary defect to issue #51
 alongside PR #56

CodeRabbit review: the CHANGELOG attributes the defect to #51, so the
active-interpolation THD comment now names both the reporting issue and
the fixing PR instead of the PR alone.
---
 internal/engine/quality_regression_test.go | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/internal/engine/quality_regression_test.go b/internal/engine/quality_regression_test.go
index ff0bb11..b3082c9 100644
--- a/internal/engine/quality_regression_test.go
+++ b/internal/engine/quality_regression_test.go
@@ -165,9 +165,9 @@ func TestQualityRegression_THD(t *testing.T) {
 // Active-interpolation THD regression (#54): unlike the exact-rational
 // ratios above, these ratios consult the sub-phase Catmull-Rom
 // interpolation banks, so this pins the whole pipeline against the
-// phase-boundary indexing defect fixed in #56 (stage-level guard:
-// phase_wrap_measure_test.go). Thresholds are measured values plus margin;
-// see constants.
+// phase-boundary indexing defect from #51, fixed in PR #56 (stage-level
+// guard: phase_wrap_measure_test.go). Thresholds are measured values plus
+// margin; see constants.
 func TestQualityRegression_THD_ActiveInterpolation(t *testing.T) {
 	tests := []struct {
 		inputRate, outputRate float64