Merge branch 'sbc-optimisation'

COMMIT_MSG.txt

@@ -1,83 +0,0 @@
-fix: 13 bugs from systematic codebase review (fix25)
-
-=== SIGNAL PATH (fixes audible/measurable on HW) ===
-
-[CRITICAL] stereo: fix 38kHz subcarrier 1-sample phase offset
-  Encode() called Sample() before capturing Phase(), so the 38kHz
-  subcarrier used phase_{n+1} while pilot used phase_n — a constant
-  60° phase error at 228kHz, degrading stereo separation to ~50%.
-  Now captures phase BEFORE tick; stores lastPhase for coherent
-  RDS carrier derivation.
-
-[HIGH] rds: phase-lock 57kHz carrier to pilot via StereoEncoder
-  RDS encoder had its own free-running 57kHz oscillator instead of
-  deriving from 3×pilot. Two independent float64 oscillators drift
-  apart over hours. Added NextSampleWithCarrier(carrier float64),
-  generator now passes stereoEncoder.RDSCarrier() (sin(3·pilotPhase))
-  so all subcarriers share one phase reference. NextSample() remains
-  as backward-compat fallback with internal carrier.
-
-[MEDIUM] dsp/fmmod: fix phase wrapping for negative phase
-  Wrap condition only triggered on phase > pi. Negative phase from
-  negative-biased composite could grow unbounded, losing float64
-  precision after extended runtime. Now wraps on |phase| > pi.
-
-[MEDIUM] dsp/oscillator: phase wrapping handles both directions
-  Same pattern as fmmod — only wrapped positive overflow. Now wraps
-  on phase >= 1 OR phase < 0 for defensive robustness.
-
-=== DAUERLAUF / SBC STABILITY ===
-
-[HIGH] offline/generator: eliminate per-chunk allocation (912KB @ 2.28MHz)
-  GenerateFrame() allocated a new []IQSample every call. At Pluto
-  rate (2.28MHz, 50ms chunks) that's 114k*8=912KB * 20/sec = 18MB/s
-  garbage, causing GC pauses and hardware underruns on Raspi.
-  Now pre-allocates frameBuf, reuses across calls. Grow-only policy,
-  never shrinks. Safe because driver.Write() is blocking.
-
-[MEDIUM] output/file: batch-write entire frame in one syscall
-  Was writing 8 bytes per sample = 114k syscalls per chunk on Pluto.
-  Now serializes into a reusable byte buffer, single file.Write().
-  Orders of magnitude faster on Raspi SD-card I/O.
-
-[MEDIUM] app/engine: add error backoff in TX loop
-  run() tight-looped at 100% CPU on persistent driver errors
-  (hardware disconnect, USB reset). Now backs off by chunkDuration
-  per error using select with ctx.Done() for clean cancellation.
-
-[MEDIUM] app/engine: replace time.Sleep shutdown with sync.WaitGroup
-  Stop() used time.Sleep(2*chunkDuration) hoping run() would exit.
-  Race condition if hardware Write() stalls. Now uses wg.Wait() for
-  deterministic goroutine join before Flush/Stop.
-
-=== CORRECTNESS / HYGIENE ===
-
-[LOW] rds/normalize: filter to ASCII, prevent mid-rune truncation
-  normalizePS truncated at 8 bytes, not 8 characters. UTF-8 input
-  (e.g. umlauts) could split mid-rune producing corrupt RDS bitstream.
-  Now replaces non-ASCII with space before truncation.
-
-[LOW] offline/generator: increment frame sequence counter
-  Sequence was hardcoded to 1 on every frame, useless for debugging.
-  Now increments per GenerateFrame() call.
-
-[LOW] control: document that config PATCH doesn't reach running engine
-  Added TODO noting that POST /config updates server's copy only;
-  running Engine/Generator holds its own snapshot.
-
-[COSMETIC] dsp/preemphasis: remove dead fields y1, a1
-  PreEmphasis.y1 and .a1 were never read in Process(). Removed from
-  struct, constructor, and Reset(). Fixed misleading doc comment on
-  PreEmphasizedSource (claims audio-rate, actually composite-rate).
-
-Files changed:
-  internal/stereo/encoder.go      - phase coherence fix
-  internal/rds/encoder.go         - pilot-locked carrier API
-  internal/rds/normalize.go       - ASCII safety
-  internal/offline/generator.go   - buffer reuse + sequence + doc
-  internal/dsp/fmmod.go           - bidirectional phase wrap
-  internal/dsp/oscillator.go      - bidirectional phase wrap
-  internal/dsp/preemphasis.go     - dead field removal
-  internal/output/file.go         - batch write
-  internal/app/engine.go          - WaitGroup + backoff
-  internal/control/control.go     - TODO config propagation

internal/app/engine.go

@@ -3,11 +3,13 @@ package app
 import (
 	"context"
 	"fmt"
+	"log"
 	"sync"
 	"sync/atomic"
 	"time"
 
 	cfgpkg "github.com/jan/fm-rds-tx/internal/config"
+	"github.com/jan/fm-rds-tx/internal/dsp"
 	offpkg "github.com/jan/fm-rds-tx/internal/offline"
 	"github.com/jan/fm-rds-tx/internal/platform"
 )
@@ -51,6 +53,7 @@ type Engine struct {
 	cfg           cfgpkg.Config
 	driver        platform.SoapyDriver
 	generator     *offpkg.Generator
+	upsampler     *dsp.FMUpsampler // nil = same-rate, non-nil = split-rate
 	chunkDuration time.Duration
 	deviceRate    float64
 
@@ -67,18 +70,41 @@ type Engine struct {
 }
 
 func NewEngine(cfg cfgpkg.Config, driver platform.SoapyDriver) *Engine {
-	// Run entire DSP chain at device rate. RDS encoder resamples its
-	// PiFmRds waveform internally. No phase upsampling needed.
 	deviceRate := cfg.EffectiveDeviceRate()
-	if deviceRate > 0 {
-		cfg.FM.CompositeRateHz = int(deviceRate)
+	compositeRate := float64(cfg.FM.CompositeRateHz)
+	if compositeRate <= 0 {
+		compositeRate = 228000
+	}
+
+	var upsampler *dsp.FMUpsampler
+
+	if deviceRate > compositeRate*1.001 {
+		// Split-rate: DSP chain runs at compositeRate (typ. 228 kHz),
+		// FMUpsampler handles FM modulation + interpolation to deviceRate.
+		// This halves CPU load compared to running all DSP at deviceRate.
+		cfg.FM.FMModulationEnabled = false
+		maxDev := cfg.FM.MaxDeviationHz
+		if maxDev <= 0 {
+			maxDev = 75000
+		}
+		upsampler = dsp.NewFMUpsampler(compositeRate, deviceRate, maxDev)
+		log.Printf("engine: split-rate mode — DSP@%.0fHz → upsample@%.0fHz (ratio %.2f)",
+			compositeRate, deviceRate, deviceRate/compositeRate)
+	} else {
+		// Same-rate: entire DSP chain runs at deviceRate.
+		// Used when deviceRate ≈ compositeRate (e.g. LimeSDR at 228 kHz).
+		if deviceRate > 0 {
+			cfg.FM.CompositeRateHz = int(deviceRate)
+		}
+		cfg.FM.FMModulationEnabled = true
+		log.Printf("engine: same-rate mode — DSP@%dHz", cfg.FM.CompositeRateHz)
 	}
-	cfg.FM.FMModulationEnabled = true
 
 	return &Engine{
 		cfg:           cfg,
 		driver:        driver,
 		generator:     offpkg.NewGenerator(cfg),
+		upsampler:     upsampler,
 		chunkDuration: 50 * time.Millisecond,
 		deviceRate:    deviceRate,
 		state:         EngineIdle,
@@ -167,6 +193,9 @@ func (e *Engine) run(ctx context.Context) {
 			return
 		}
 		frame := e.generator.GenerateFrame(e.chunkDuration)
+		if e.upsampler != nil {
+			frame = e.upsampler.Process(frame)
+		}
 		n, err := e.driver.Write(ctx, frame)
 		if err != nil {
 			if ctx.Err() != nil { return }

internal/app/engine_test.go

@@ -79,3 +79,55 @@ func TestEngineDriverStats(t *testing.T) {
 		t.Fatal("expected driver to have written frames")
 	}
 }
+
+func TestEngineSplitRate(t *testing.T) {
+	// Simulate Pluto-like config: composite at 228kHz, device at 2.28MHz
+	cfg := cfgpkg.Default()
+	cfg.Backend.DeviceSampleRateHz = 2280000 // 10:1 ratio
+
+	driver := platform.NewSimulatedDriver(nil)
+	eng := NewEngine(cfg, driver)
+	eng.SetChunkDuration(10 * time.Millisecond)
+
+	// Verify split-rate path was chosen
+	if eng.upsampler == nil {
+		t.Fatal("expected split-rate mode (upsampler != nil)")
+	}
+
+	ctx := context.Background()
+	if err := eng.Start(ctx); err != nil {
+		t.Fatalf("start: %v", err)
+	}
+
+	time.Sleep(200 * time.Millisecond)
+
+	stats := eng.Stats()
+	if stats.ChunksProduced < 5 {
+		t.Fatalf("expected at least 5 chunks, got %d", stats.ChunksProduced)
+	}
+	// With 10:1 upsampling and 10ms chunks at 228kHz source:
+	// 2280 src samples → 22800 output samples per chunk
+	// 5 chunks minimum → at least 114000 samples
+	if stats.TotalSamples < 50000 {
+		t.Fatalf("expected at least 50000 upsampled samples, got %d", stats.TotalSamples)
+	}
+
+	if err := eng.Stop(ctx); err != nil {
+		t.Fatalf("stop: %v", err)
+	}
+	if stats.Underruns > 0 {
+		t.Fatalf("unexpected underruns: %d", stats.Underruns)
+	}
+}
+
+func TestEngineSameRate(t *testing.T) {
+	// Default config: deviceSampleRateHz=0, compositeRateHz=228000
+	// EffectiveDeviceRate = 228000 → same-rate mode, no upsampler
+	cfg := cfgpkg.Default()
+	driver := platform.NewSimulatedDriver(nil)
+	eng := NewEngine(cfg, driver)
+
+	if eng.upsampler != nil {
+		t.Fatal("expected same-rate mode (upsampler == nil)")
+	}
+}

internal/dsp/fmupsample.go

@@ -0,0 +1,189 @@
+package dsp
+
+import (
+	"math"
+
+	"github.com/jan/fm-rds-tx/internal/output"
+)
+
+// FMUpsampler converts a composite baseband signal at a low source rate into
+// FM-modulated IQ samples at a higher device rate, via phase-domain
+// interpolation.
+//
+// Architecture: accumulate FM phase at source rate (cheap, few trig ops),
+// then linearly interpolate the phase to device rate and emit sin/cos.
+// This is mathematically equivalent to running the full FMModulator at device
+// rate, but needs trig only at the output rate — saving all the DSP that
+// would otherwise run at the higher rate (stereo encode, RDS, limiter etc.).
+//
+// Cross-chunk boundary: the upsampler carries over prevPhase and srcPos
+// between calls. The interpolation coordinate system places prevPhase at
+// virtual index 0 and srcPhases[0..N-1] at indices 1..N. This guarantees
+// smooth phase transitions at every chunk boundary with zero discontinuity.
+//
+// Zero-allocation in steady state: all buffers are pre-allocated on first
+// call and reused. The returned CompositeFrame is an internal buffer —
+// valid only until the next Process() call.
+type FMUpsampler struct {
+	srcRate      float64 // composite rate (e.g. 228000)
+	dstRate      float64 // device rate (e.g. 2280000)
+	maxDeviation float64 // peak FM deviation in Hz (e.g. 75000)
+	step         float64 // source-samples per output-sample = srcRate/dstRate
+
+	// Persistent state across Process() calls
+	phase     float64 // accumulated FM phase in radians, continuous across chunks
+	prevPhase float64 // phase at end of previous chunk (virtual index 0)
+	srcPos    float64 // fractional source position carry-over into next chunk
+	seeded    bool    // true after first Process() call
+
+	// Pre-allocated buffers — grown once, never shrunk
+	srcPhases []float64
+	outBuf    []output.IQSample
+	outFrame  output.CompositeFrame
+}
+
+// NewFMUpsampler creates a phase-domain upsampler.
+//
+// srcRate:      composite DSP rate (typ. 228000 Hz)
+// dstRate:      device output rate (typ. 2280000 Hz)
+// maxDeviation: FM peak deviation (typ. 75000 Hz)
+func NewFMUpsampler(srcRate, dstRate, maxDeviation float64) *FMUpsampler {
+	return &FMUpsampler{
+		srcRate:      srcRate,
+		dstRate:      dstRate,
+		maxDeviation: maxDeviation,
+		step:         srcRate / dstRate,
+	}
+}
+
+// Process takes a CompositeFrame where Samples[i].I contains the composite
+// baseband value (FM modulation must be OFF in the generator; Q is ignored).
+// Returns an FM-modulated IQ frame at dstRate.
+//
+// The returned frame is an internal buffer — valid until the next Process()
+// call. The caller must consume or copy the data before calling again.
+func (u *FMUpsampler) Process(frame *output.CompositeFrame) *output.CompositeFrame {
+	if frame == nil || len(frame.Samples) == 0 {
+		return frame
+	}
+
+	srcLen := len(frame.Samples)
+
+	// --- Phase accumulation at source rate ---
+
+	// Grow srcPhases buffer if needed
+	if cap(u.srcPhases) < srcLen {
+		u.srcPhases = make([]float64, srcLen)
+	}
+	srcPhases := u.srcPhases[:srcLen]
+
+	phaseInc := 2 * math.Pi * u.maxDeviation / u.srcRate
+	for i, s := range frame.Samples {
+		u.phase += float64(s.I) * phaseInc
+		srcPhases[i] = u.phase
+	}
+
+	// Phase wrapping — symmetric, shift prevPhase in lockstep
+	if u.phase > math.Pi || u.phase < -math.Pi {
+		offset := 2 * math.Pi * math.Floor((u.phase+math.Pi)/(2*math.Pi))
+		u.phase -= offset
+		for i := range srcPhases {
+			srcPhases[i] -= offset
+		}
+		if u.seeded {
+			u.prevPhase -= offset
+		}
+	}
+
+	// Seed prevPhase on very first call
+	if !u.seeded {
+		// Extrapolate backwards from first two phases to get a virtual "previous"
+		// phase, so the first chunk's boundary interpolation is smooth.
+		if srcLen >= 2 {
+			u.prevPhase = 2*srcPhases[0] - srcPhases[1]
+		} else {
+			u.prevPhase = srcPhases[0]
+		}
+		u.srcPos = 0
+		u.seeded = true
+	}
+
+	// --- Interpolation coordinate system ---
+	//
+	// Virtual index 0 = prevPhase (end of previous chunk)
+	// Virtual index 1 = srcPhases[0]
+	// Virtual index 2 = srcPhases[1]
+	// ...
+	// Virtual index N = srcPhases[N-1]
+	//
+	// srcPos ranges from 0 (carry-over) to N (= srcLen).
+	// We generate output samples while srcPos < srcLen (virtual index srcLen).
+
+	// phaseAt returns the phase at a virtual index.
+	phaseAt := func(vi int) float64 {
+		if vi <= 0 {
+			return u.prevPhase
+		}
+		if vi > srcLen {
+			return srcPhases[srcLen-1]
+		}
+		return srcPhases[vi-1]
+	}
+
+	// Calculate output count: from srcPos to srcLen, stepping by u.step.
+	// +1 for safety margin; we clamp below.
+	maxOut := int(math.Ceil(float64(srcLen)-u.srcPos)/u.step) + 1
+	if maxOut < 0 {
+		maxOut = 0
+	}
+
+	// Grow output buffer if needed
+	if cap(u.outBuf) < maxOut {
+		u.outBuf = make([]output.IQSample, maxOut)
+	}
+	out := u.outBuf[:maxOut]
+
+	// --- Generate output samples ---
+	pos := u.srcPos
+	n := 0
+	for pos < float64(srcLen) && n < maxOut {
+		vi := int(pos)     // virtual index (integer part)
+		frac := pos - float64(vi)
+
+		pA := phaseAt(vi)
+		pB := phaseAt(vi + 1)
+		p := pA + frac*(pB-pA)
+
+		out[n] = output.IQSample{
+			I: float32(math.Cos(p)),
+			Q: float32(math.Sin(p)),
+		}
+		n++
+		pos += u.step
+	}
+
+	// Carry state for next chunk
+	u.prevPhase = srcPhases[srcLen-1]
+	u.srcPos = pos - float64(srcLen)
+
+	// Package output
+	u.outFrame.Samples = out[:n]
+	u.outFrame.SampleRateHz = u.dstRate
+	u.outFrame.Timestamp = frame.Timestamp
+	u.outFrame.Sequence = frame.Sequence
+
+	return &u.outFrame
+}
+
+// Reset clears all accumulated state, as if freshly constructed.
+func (u *FMUpsampler) Reset() {
+	u.phase = 0
+	u.prevPhase = 0
+	u.srcPos = 0
+	u.seeded = false
+}
+
+// Stats returns internal state for diagnostics/testing.
+func (u *FMUpsampler) Stats() (phase, prevPhase, srcPos float64) {
+	return u.phase, u.prevPhase, u.srcPos
+}

internal/dsp/fmupsample_test.go

@@ -0,0 +1,617 @@
+package dsp
+
+import (
+	"math"
+	"testing"
+
+	"github.com/jan/fm-rds-tx/internal/output"
+)
+
+// makeCompositeFrame builds a CompositeFrame with composite values in I, Q=0.
+func makeCompositeFrame(values []float64, rate float64) *output.CompositeFrame {
+	samples := make([]output.IQSample, len(values))
+	for i, v := range values {
+		samples[i] = output.IQSample{I: float32(v), Q: 0}
+	}
+	return &output.CompositeFrame{
+		Samples:      samples,
+		SampleRateHz: rate,
+		Sequence:     1,
+	}
+}
+
+// iqMagnitude returns sqrt(I²+Q²) for an IQ sample.
+func iqMagnitude(s output.IQSample) float64 {
+	return math.Sqrt(float64(s.I)*float64(s.I) + float64(s.Q)*float64(s.Q))
+}
+
+// iqPhase returns atan2(Q, I) for an IQ sample.
+func iqPhase(s output.IQSample) float64 {
+	return math.Atan2(float64(s.Q), float64(s.I))
+}
+
+// phaseDiff returns the shortest signed angular difference between two angles.
+func phaseDiff(a, b float64) float64 {
+	d := b - a
+	for d > math.Pi {
+		d -= 2 * math.Pi
+	}
+	for d < -math.Pi {
+		d += 2 * math.Pi
+	}
+	return d
+}
+
+// --- Test: zero composite produces no frequency offset ---
+
+func TestFMUpsampler_ZeroComposite(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+	n := 11400 // 50ms at 228kHz
+	vals := make([]float64, n)
+	// All zeros — FM deviation = 0, carrier should be constant phase.
+	frame := makeCompositeFrame(vals, 228000)
+	out := u.Process(frame)
+
+	if len(out.Samples) == 0 {
+		t.Fatal("no output samples")
+	}
+	if out.SampleRateHz != 2280000 {
+		t.Fatalf("expected rate 2280000, got %f", out.SampleRateHz)
+	}
+
+	// All output IQ should have magnitude ≈ 1.0 and constant phase (≈ 0)
+	for i, s := range out.Samples {
+		mag := iqMagnitude(s)
+		if math.Abs(mag-1.0) > 0.001 {
+			t.Fatalf("sample %d: magnitude %.6f, expected ~1.0", i, mag)
+		}
+	}
+
+	// Phase should not advance (zero deviation)
+	firstPhase := iqPhase(out.Samples[0])
+	for i := 1; i < len(out.Samples); i++ {
+		p := iqPhase(out.Samples[i])
+		if math.Abs(phaseDiff(firstPhase, p)) > 0.01 {
+			t.Fatalf("sample %d: phase drifted to %.4f (first was %.4f)", i, p, firstPhase)
+		}
+	}
+}
+
+// --- Test: DC composite produces constant frequency ---
+
+func TestFMUpsampler_DCComposite(t *testing.T) {
+	srcRate := 228000.0
+	dstRate := 2280000.0
+	maxDev := 75000.0
+
+	u := NewFMUpsampler(srcRate, dstRate, maxDev)
+
+	dc := 0.5 // half deviation = +37.5 kHz
+	n := 11400
+	vals := make([]float64, n)
+	for i := range vals {
+		vals[i] = dc
+	}
+	frame := makeCompositeFrame(vals, srcRate)
+	out := u.Process(frame)
+
+	if len(out.Samples) < 100 {
+		t.Fatalf("too few output samples: %d", len(out.Samples))
+	}
+
+	// Expected frequency: dc * maxDev = 37500 Hz
+	// Expected phase step per output sample: 2π * 37500 / 2280000
+	expectedStep := 2 * math.Pi * dc * maxDev / dstRate
+
+	// Measure average phase step over a stable region (skip first 100)
+	var totalDiff float64
+	count := 0
+	for i := 101; i < len(out.Samples) && i < 10000; i++ {
+		d := phaseDiff(iqPhase(out.Samples[i-1]), iqPhase(out.Samples[i]))
+		totalDiff += d
+		count++
+	}
+	avgStep := totalDiff / float64(count)
+
+	if math.Abs(avgStep-expectedStep) > 0.001 {
+		t.Fatalf("average phase step = %.6f, expected %.6f (error %.6f)",
+			avgStep, expectedStep, avgStep-expectedStep)
+	}
+}
+
+// --- Test: IQ magnitude is always ≈ 1.0 ---
+
+func TestFMUpsampler_UnitMagnitude(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	// Varying composite: a 1 kHz tone at full deviation
+	n := 11400
+	vals := make([]float64, n)
+	for i := range vals {
+		vals[i] = math.Sin(2 * math.Pi * 1000 * float64(i) / 228000)
+	}
+	frame := makeCompositeFrame(vals, 228000)
+	out := u.Process(frame)
+
+	for i, s := range out.Samples {
+		mag := iqMagnitude(s)
+		if math.Abs(mag-1.0) > 0.002 {
+			t.Fatalf("sample %d: magnitude %.6f, expected ~1.0", i, mag)
+		}
+	}
+}
+
+// --- Test: output sample count for exact integer ratio ---
+
+func TestFMUpsampler_OutputCountExactRatio(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	srcLen := 11400 // 50ms
+	vals := make([]float64, srcLen)
+	frame := makeCompositeFrame(vals, 228000)
+
+	// With exact 10:1 ratio, each chunk should produce exactly srcLen * 10 samples.
+	// Small tolerance for boundary effects on first chunk.
+	for chunk := 0; chunk < 10; chunk++ {
+		out := u.Process(frame)
+		expected := srcLen * 10
+		if out == nil || len(out.Samples) < expected-1 || len(out.Samples) > expected+1 {
+			t.Fatalf("chunk %d: got %d output samples, expected ~%d",
+				chunk, len(out.Samples), expected)
+		}
+	}
+}
+
+// --- Test: output sample count for non-integer ratio ---
+
+func TestFMUpsampler_OutputCountNonIntegerRatio(t *testing.T) {
+	// PlutoSDR minimum: 228000 → 2084000, ratio ≈ 9.14
+	u := NewFMUpsampler(228000, 2084000, 75000)
+
+	srcLen := 11400
+	vals := make([]float64, srcLen)
+	frame := makeCompositeFrame(vals, 228000)
+
+	var totalOut int
+	for chunk := 0; chunk < 20; chunk++ {
+		out := u.Process(frame)
+		totalOut += len(out.Samples)
+	}
+
+	// Over 20 chunks, total output should be close to 20 * srcLen * ratio
+	expectedTotal := int(20 * float64(srcLen) * 2084000 / 228000)
+	drift := math.Abs(float64(totalOut-expectedTotal)) / float64(expectedTotal)
+	if drift > 0.001 {
+		t.Fatalf("total output %d, expected ~%d (drift %.4f%%)", totalOut, expectedTotal, drift*100)
+	}
+}
+
+// --- CRITICAL TEST: phase continuity across chunk boundaries ---
+
+func TestFMUpsampler_PhaseContinuityAcrossChunks(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	srcLen := 11400
+	dc := 0.3 // constant deviation
+	vals := make([]float64, srcLen)
+	for i := range vals {
+		vals[i] = dc
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	var prevLastPhase float64
+	for chunk := 0; chunk < 20; chunk++ {
+		out := u.Process(frame)
+		n := len(out.Samples)
+		if n == 0 {
+			t.Fatalf("chunk %d: no output", chunk)
+		}
+
+		firstPhase := iqPhase(out.Samples[0])
+
+		// Check continuity from previous chunk
+		if chunk > 0 {
+			jump := math.Abs(phaseDiff(prevLastPhase, firstPhase))
+			// For DC composite at dstRate, the expected step between
+			// consecutive output samples is 2π * dc * maxDev / dstRate ≈ 0.062 rad.
+			// The boundary jump should be close to this, not larger.
+			expectedStep := 2 * math.Pi * dc * 75000 / 2280000
+			if jump > expectedStep*3 {
+				t.Fatalf("chunk %d: boundary phase jump %.6f rad (expected ~%.6f)",
+					chunk, jump, expectedStep)
+			}
+		}
+
+		// Also check that phase increments WITHIN the chunk are smooth
+		maxJump := 0.0
+		for i := 1; i < n; i++ {
+			d := math.Abs(phaseDiff(iqPhase(out.Samples[i-1]), iqPhase(out.Samples[i])))
+			if d > maxJump {
+				maxJump = d
+			}
+		}
+		expectedIntra := 2 * math.Pi * dc * 75000 / 2280000
+		if maxJump > expectedIntra*2 {
+			t.Fatalf("chunk %d: intra-chunk max phase jump %.6f (expected ~%.6f)",
+				chunk, maxJump, expectedIntra)
+		}
+
+		prevLastPhase = iqPhase(out.Samples[n-1])
+	}
+}
+
+// --- Test: non-integer ratio boundary continuity ---
+
+func TestFMUpsampler_BoundaryContinuityNonIntegerRatio(t *testing.T) {
+	// 228000 → 2084000, ratio ≈ 9.14 — carry-over is non-zero every chunk
+	u := NewFMUpsampler(228000, 2084000, 75000)
+
+	srcLen := 11400
+	dc := 0.5
+	vals := make([]float64, srcLen)
+	for i := range vals {
+		vals[i] = dc
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	expectedStep := 2 * math.Pi * dc * 75000 / 2084000
+	var prevLastPhase float64
+
+	for chunk := 0; chunk < 50; chunk++ {
+		out := u.Process(frame)
+		n := len(out.Samples)
+		if n == 0 {
+			t.Fatalf("chunk %d: no output", chunk)
+		}
+
+		if chunk > 0 {
+			jump := math.Abs(phaseDiff(prevLastPhase, iqPhase(out.Samples[0])))
+			if jump > expectedStep*3 {
+				t.Fatalf("chunk %d: boundary jump %.6f rad (expected ~%.6f)",
+					chunk, jump, expectedStep)
+			}
+		}
+
+		prevLastPhase = iqPhase(out.Samples[n-1])
+	}
+}
+
+// --- Test: phase wrapping doesn't introduce discontinuity ---
+
+func TestFMUpsampler_PhaseWrapping(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	// Full deviation for many chunks — phase will wrap many times
+	srcLen := 11400
+	vals := make([]float64, srcLen)
+	for i := range vals {
+		vals[i] = 1.0 // full positive deviation
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	expectedStep := 2 * math.Pi * 1.0 * 75000 / 2280000
+	var prevLastPhase float64
+
+	for chunk := 0; chunk < 100; chunk++ {
+		out := u.Process(frame)
+		n := len(out.Samples)
+
+		if chunk > 0 {
+			jump := math.Abs(phaseDiff(prevLastPhase, iqPhase(out.Samples[0])))
+			if jump > expectedStep*3 {
+				t.Fatalf("chunk %d: post-wrap boundary jump %.6f (limit %.6f)",
+					chunk, jump, expectedStep*3)
+			}
+		}
+
+		// Spot-check magnitudes (should still be 1.0 after wrapping)
+		for i := 0; i < n; i += n / 10 {
+			mag := iqMagnitude(out.Samples[i])
+			if math.Abs(mag-1.0) > 0.002 {
+				t.Fatalf("chunk %d sample %d: magnitude %.6f after wrapping", chunk, i, mag)
+			}
+		}
+
+		prevLastPhase = iqPhase(out.Samples[n-1])
+	}
+}
+
+// --- Test: negative composite (phase wraps negative direction) ---
+
+func TestFMUpsampler_NegativeDeviation(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	srcLen := 11400
+	vals := make([]float64, srcLen)
+	for i := range vals {
+		vals[i] = -0.8
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	expectedStep := 2 * math.Pi * 0.8 * 75000 / 2280000 // magnitude
+	var prevLastPhase float64
+
+	for chunk := 0; chunk < 100; chunk++ {
+		out := u.Process(frame)
+		n := len(out.Samples)
+
+		if chunk > 0 {
+			jump := math.Abs(phaseDiff(prevLastPhase, iqPhase(out.Samples[0])))
+			if jump > expectedStep*3 {
+				t.Fatalf("chunk %d: negative-dev boundary jump %.6f", chunk, jump)
+			}
+		}
+		prevLastPhase = iqPhase(out.Samples[n-1])
+	}
+}
+
+// --- Test: equivalence with direct FMModulator (frequency comparison) ---
+
+func TestFMUpsampler_EquivalenceWithFMModulator(t *testing.T) {
+	srcRate := 228000.0
+	maxDev := 75000.0
+
+	// Generate a 1kHz tone composite signal
+	srcLen := 11400
+	composite := make([]float64, srcLen)
+	for i := range composite {
+		composite[i] = 0.6 * math.Sin(2*math.Pi*1000*float64(i)/srcRate)
+	}
+
+	// Path A: FMUpsampler at 1:1 ratio (no upsampling, just FM modulation)
+	up := NewFMUpsampler(srcRate, srcRate, maxDev)
+	frame := makeCompositeFrame(composite, srcRate)
+	outA := up.Process(frame)
+
+	// Path B: Direct FMModulator
+	mod := NewFMModulator(srcRate)
+	mod.MaxDeviation = maxDev
+	outB := make([]output.IQSample, srcLen)
+	for i, c := range composite {
+		oi, oq := mod.Modulate(c)
+		outB[i] = output.IQSample{I: float32(oi), Q: float32(oq)}
+	}
+
+	// The upsampler has a 1-sample offset due to the virtual prevPhase index.
+	// Instead of comparing absolute phase, compare instantaneous frequency
+	// (phase step per sample), which is invariant to time shifts.
+
+	freqA := make([]float64, 0, len(outA.Samples)-1)
+	for i := 1; i < len(outA.Samples); i++ {
+		freqA = append(freqA, phaseDiff(iqPhase(outA.Samples[i-1]), iqPhase(outA.Samples[i])))
+	}
+
+	freqB := make([]float64, 0, len(outB)-1)
+	for i := 1; i < len(outB); i++ {
+		freqB = append(freqB, phaseDiff(iqPhase(outB[i-1]), iqPhase(outB[i])))
+	}
+
+	// Compare frequency profiles (skip edges, allow 1-sample alignment offset)
+	minLen := len(freqA)
+	if len(freqB) < minLen {
+		minLen = len(freqB)
+	}
+
+	maxFreqDiff := 0.0
+	for i := 20; i < minLen-20; i++ {
+		// Try both aligned and 1-sample-shifted
+		d0 := math.Abs(freqA[i] - freqB[i])
+		d1 := math.Abs(freqA[i] - freqB[i-1])
+		d := math.Min(d0, d1)
+		if d > maxFreqDiff {
+			maxFreqDiff = d
+		}
+	}
+
+	// Instantaneous frequency should match within 0.02 rad/sample
+	if maxFreqDiff > 0.02 {
+		t.Fatalf("max instantaneous frequency difference: %.6f rad/sample (too large)", maxFreqDiff)
+	}
+
+	// All magnitudes should be ~1.0
+	maxMagDiff := 0.0
+	for i := range outA.Samples {
+		mag := iqMagnitude(outA.Samples[i])
+		d := math.Abs(mag - 1.0)
+		if d > maxMagDiff {
+			maxMagDiff = d
+		}
+	}
+	if maxMagDiff > 0.01 {
+		t.Fatalf("max magnitude deviation: %.6f", maxMagDiff)
+	}
+
+	t.Logf("equivalence: maxFreqDiff=%.6f rad/sample, maxMagDiff=%.6f", maxFreqDiff, maxMagDiff)
+}
+
+// --- Test: buffer reuse (no allocation after first call) ---
+
+func TestFMUpsampler_BufferReuse(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	srcLen := 11400
+	vals := make([]float64, srcLen)
+	frame := makeCompositeFrame(vals, 228000)
+
+	// First call allocates
+	out1 := u.Process(frame)
+	ptr1 := &out1.Samples[0]
+
+	// Second call should reuse the same buffer
+	out2 := u.Process(frame)
+	ptr2 := &out2.Samples[0]
+
+	if ptr1 != ptr2 {
+		t.Fatal("buffer was reallocated on second call — should reuse")
+	}
+}
+
+// --- Test: Reset clears state properly ---
+
+func TestFMUpsampler_Reset(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	vals := make([]float64, 11400)
+	for i := range vals {
+		vals[i] = 0.5
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	// Run a few chunks to accumulate state
+	for i := 0; i < 5; i++ {
+		u.Process(frame)
+	}
+
+	phase1, prev1, pos1 := u.Stats()
+	if phase1 == 0 && prev1 == 0 && pos1 == 0 {
+		t.Fatal("state should be non-zero after processing")
+	}
+
+	u.Reset()
+	phase2, prev2, pos2 := u.Stats()
+	if phase2 != 0 || prev2 != 0 || pos2 != 0 {
+		t.Fatalf("state not zeroed after Reset: phase=%f prev=%f pos=%f", phase2, prev2, pos2)
+	}
+}
+
+// --- Test: tiny chunks (edge case: 1 source sample) ---
+
+func TestFMUpsampler_TinyChunk(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	vals := []float64{0.5}
+	frame := makeCompositeFrame(vals, 228000)
+
+	for chunk := 0; chunk < 20; chunk++ {
+		out := u.Process(frame)
+		if len(out.Samples) == 0 {
+			t.Fatalf("chunk %d: no output from single-sample input", chunk)
+		}
+		for i, s := range out.Samples {
+			mag := iqMagnitude(s)
+			if math.Abs(mag-1.0) > 0.01 {
+				t.Fatalf("chunk %d sample %d: magnitude %.4f", chunk, i, mag)
+			}
+		}
+	}
+}
+
+// --- Test: alternating sign composite (stress boundary interpolation) ---
+
+func TestFMUpsampler_AlternatingComposite(t *testing.T) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	srcLen := 11400
+	vals := make([]float64, srcLen)
+	for i := range vals {
+		if i%2 == 0 {
+			vals[i] = 0.8
+		} else {
+			vals[i] = -0.8
+		}
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	var prevLastPhase float64
+	for chunk := 0; chunk < 20; chunk++ {
+		out := u.Process(frame)
+		n := len(out.Samples)
+
+		if chunk > 0 {
+			jump := math.Abs(phaseDiff(prevLastPhase, iqPhase(out.Samples[0])))
+			// Alternating composite: the phase meanders. Boundary jump
+			// should still be bounded by the maximum single-sample phase change.
+			maxSingleStep := 2 * math.Pi * 0.8 * 75000 / 228000.0 // at source rate
+			maxOutputJump := maxSingleStep * 0.1                    // step = srcRate/dstRate
+			if jump > maxOutputJump*3 {
+				t.Fatalf("chunk %d: alternating boundary jump %.4f (limit %.4f)",
+					chunk, jump, maxOutputJump*3)
+			}
+		}
+
+		// Check all magnitudes
+		for i, s := range out.Samples {
+			mag := iqMagnitude(s)
+			if math.Abs(mag-1.0) > 0.01 {
+				t.Fatalf("chunk %d sample %d: magnitude %.4f", chunk, i, mag)
+			}
+		}
+
+		prevLastPhase = iqPhase(out.Samples[n-1])
+	}
+}
+
+// --- Test: long-running stability (simulates 1 hour) ---
+
+func TestFMUpsampler_LongRun(t *testing.T) {
+	if testing.Short() {
+		t.Skip("skipping long-run test in short mode")
+	}
+
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	srcLen := 11400
+	vals := make([]float64, srcLen)
+	for i := range vals {
+		vals[i] = 0.3 * math.Sin(2*math.Pi*1000*float64(i)/228000)
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	// 20 chunks/sec × 3600 sec = 72000 chunks per hour.
+	// We test 72000 chunks (simulating 1 hour).
+	chunks := 72000
+	var prevLastPhase float64
+
+	for chunk := 0; chunk < chunks; chunk++ {
+		out := u.Process(frame)
+		n := len(out.Samples)
+
+		if n == 0 {
+			t.Fatalf("chunk %d: no output", chunk)
+		}
+
+		// Check boundary every 1000 chunks
+		if chunk > 0 && chunk%1000 == 0 {
+			jump := math.Abs(phaseDiff(prevLastPhase, iqPhase(out.Samples[0])))
+			if jump > 0.5 {
+				t.Fatalf("chunk %d: boundary jump %.4f rad (phase drift?)", chunk, jump)
+			}
+
+			// Check magnitude of first sample
+			mag := iqMagnitude(out.Samples[0])
+			if math.Abs(mag-1.0) > 0.01 {
+				t.Fatalf("chunk %d: magnitude %.4f after long run", chunk, mag)
+			}
+		}
+
+		prevLastPhase = iqPhase(out.Samples[n-1])
+	}
+
+	// Check phase is still bounded
+	phase, _, _ := u.Stats()
+	if math.Abs(phase) > 2*math.Pi {
+		t.Fatalf("phase unbounded after %d chunks: %.2f", chunks, phase)
+	}
+}
+
+// --- Benchmark ---
+
+func BenchmarkFMUpsampler_Process(b *testing.B) {
+	u := NewFMUpsampler(228000, 2280000, 75000)
+
+	srcLen := 11400
+	vals := make([]float64, srcLen)
+	for i := range vals {
+		vals[i] = 0.5 * math.Sin(2*math.Pi*1000*float64(i)/228000)
+	}
+	frame := makeCompositeFrame(vals, 228000)
+
+	// Warm up (trigger buffer allocation)
+	u.Process(frame)
+
+	b.ResetTimer()
+	b.ReportAllocs()
+	for i := 0; i < b.N; i++ {
+		u.Process(frame)
+	}
+}