fm-rds-tx/internal/dsp/compositeclipper.go

// Package dsp — CompositeClipper implements a broadcast-grade iterative
// composite MPX clipper with optional soft-knee clipping and look-ahead
// peak limiting.
//
// # Signal flow
//
//   audioMPX (mono + stereo sub, no pilot/RDS)
//     → Look-ahead pre-limiter (optional, reduces peak:average before clipping)
//     → N iterations of: Clip → Notch19kHz → Notch57kHz
//     → Final hard-clip safety net
//     → output
//
// Each iteration reduces spectral splatter into the pilot (19 kHz) and RDS
// (57 kHz) guard bands caused by the previous clip, while the subsequent
// clip catches the filter overshoot. After 2-3 iterations the signal
// converges: peaks stay below ceiling AND protected bands are clean.
//
// The look-ahead limiter uses a delay line so the envelope detector sees
// peaks before they reach the clipper. This allows smooth gain reduction
// instead of hard clipping, producing fewer harmonics per iteration.
//
// # Comparable to
//
// Omnia.11 composite processing, Orban Optimod 8700 composite clipper,
// and similar ITU-R SM.1268-compliant broadcast processors.
package dsp

import "math"

// CompositeClipperConfig holds the parameters for NewCompositeClipper.
type CompositeClipperConfig struct {
	// Ceiling is the peak output level. Typically 1.0 (= 100% modulation
	// for the audio portion). Pilot and RDS are added externally after.
	Ceiling float64

	// Iterations is the number of clip-filter-clip passes (1-5).
	// More iterations = cleaner guard bands, slightly more latency from
	// filter group delay. 3 is a good default; 1 is conservative, 5 is
	// aggressive (Omnia "brick wall" territory).
	Iterations int

	// SoftKnee sets the width of the soft-clip transition zone below
	// ceiling, in linear amplitude. 0 = hard clip, 0.15 = moderate,
	// 0.3 = gentle. The soft clipper uses tanh() compression above
	// (ceiling - softKnee), producing fewer harmonics per iteration
	// at the cost of slightly lower peak density.
	SoftKnee float64

	// LookaheadMs sets the look-ahead delay in milliseconds. 0 = disabled.
	// Typical values: 0.5-2.0 ms. The limiter pre-reduces peaks before
	// they hit the clipper, so the clipper does less work and generates
	// fewer harmonics. Introduces LookaheadMs of audio latency.
	LookaheadMs float64

	// SampleRate must match the composite DSP rate (typically 228000 Hz).
	SampleRate float64
}

// CompositeClipper is a stateful per-sample processor. Not thread-safe —
// call Process from the single DSP goroutine only.
type CompositeClipper struct {
	ceiling    float64
	softKnee   float64
	iterations int

	// Per-iteration filter banks — each iteration needs independent state
	// because it processes a different signal (output of previous clip).
	notch19 []*FilterChain
	notch57 []*FilterChain

	// Look-ahead limiter state
	laEnabled bool
	laDelay   []float64 // circular delay buffer
	laLen     int       // delay length in samples
	laWrite   int       // write position
	laEnv     float64   // peak envelope (tracks input peaks)
	laGain    float64   // current gain multiplier (applied to delayed output)
	laAttack  float64   // envelope/gain attack coefficient
	laRelease float64   // envelope/gain release coefficient
}

// NewCompositeClipper creates a ready-to-use composite clipper.
// Returns nil if cfg is invalid or sample rate is zero.
func NewCompositeClipper(cfg CompositeClipperConfig) *CompositeClipper {
	if cfg.SampleRate <= 0 {
		return nil
	}
	if cfg.Ceiling <= 0 {
		cfg.Ceiling = 1.0
	}
	if cfg.Iterations < 1 {
		cfg.Iterations = 1
	}
	if cfg.Iterations > 5 {
		cfg.Iterations = 5
	}
	if cfg.SoftKnee < 0 {
		cfg.SoftKnee = 0
	}

	c := &CompositeClipper{
		ceiling:    cfg.Ceiling,
		softKnee:   cfg.SoftKnee,
		iterations: cfg.Iterations,
		notch19:    make([]*FilterChain, cfg.Iterations),
		notch57:    make([]*FilterChain, cfg.Iterations),
		laGain:     1.0,
	}

	// Build per-iteration filter pairs.
	// Double-cascade notch at each frequency for deep rejection.
	// Q=15 at 19 kHz → narrow (~1.3 kHz), preserves stereo sub.
	// Q=10 at 57 kHz → slightly wider (~5.7 kHz), covers RDS bandwidth.
	for i := 0; i < cfg.Iterations; i++ {
		c.notch19[i] = &FilterChain{
			Stages: []Biquad{
				*NewNotch(19000, cfg.SampleRate, 15),
				*NewNotch(19000, cfg.SampleRate, 15),
			},
		}
		c.notch57[i] = &FilterChain{
			Stages: []Biquad{
				*NewNotch(57000, cfg.SampleRate, 10),
				*NewNotch(57000, cfg.SampleRate, 10),
			},
		}
	}

	// Look-ahead limiter
	if cfg.LookaheadMs > 0 {
		c.laEnabled = true
		c.laLen = int(math.Round(cfg.LookaheadMs * cfg.SampleRate / 1000))
		if c.laLen < 1 {
			c.laLen = 1
		}
		c.laDelay = make([]float64, c.laLen)

		// Attack: ramp down over half the look-ahead window so gain is
		// fully reduced by the time the peak exits the delay line.
		// Using half the window gives ~86% reduction at exit (2 time constants).
		attackSamples := float64(c.laLen) / 2
		if attackSamples < 1 {
			attackSamples = 1
		}
		c.laAttack = 1.0 - math.Exp(-1.0/attackSamples)

		// Release: 20ms — fast enough to recover between transients,
		// slow enough to avoid gain pumping on dense material.
		releaseSamples := 0.020 * cfg.SampleRate
		if releaseSamples < 1 {
			releaseSamples = 1
		}
		c.laRelease = 1.0 - math.Exp(-1.0/releaseSamples)
	}

	return c
}

// Process runs one composite sample through the full chain.
//
// Input: audio-only MPX (mono + stereo sub, no pilot, no RDS).
// Output: peak-limited, spectrally clean MPX. Never exceeds ceiling.
//
// The caller adds pilot and RDS to the output AFTER this call.
func (c *CompositeClipper) Process(audioMPX float64) float64 {
	x := audioMPX

	// --- Stage 1: Look-ahead pre-limiting ---
	// Reduces peaks smoothly before the clipper sees them.
	// Fewer hard-clip events → fewer harmonics → cleaner output.
	if c.laEnabled {
		x = c.processLookahead(x)
	}

	// --- Stage 2: Iterative clip → notch → notch ---
	// Each pass:
	//  - Clip catches peaks (from input or previous filter overshoot)
	//  - Notch19 removes energy at pilot frequency
	//  - Notch57 removes energy at RDS frequency
	//  - Filter ringing creates small overshoots → next pass catches them
	// After N iterations the overshoot converges to near-zero.
	for i := 0; i < c.iterations; i++ {
		x = c.clip(x)
		x = c.notch19[i].Process(x)
		x = c.notch57[i].Process(x)
	}

	// --- Stage 3: Final safety hard-clip ---
	// Catches residual overshoot from the last iteration's notch filters.
	// Even with 3+ iterations there can be sub-0.1 dB overshoot from the
	// final notch; this guarantees the output never exceeds ceiling.
	x = HardClip(x, c.ceiling)

	return x
}

// processLookahead applies look-ahead peak limiting.
//
// The input enters a delay line. The envelope detector runs on the INPUT
// (not the delayed output), so it sees peaks LookaheadMs before they exit
// the delay. Gain is reduced smoothly over the look-ahead window, so when
// the peak arrives at the output, the gain is already low enough.
func (c *CompositeClipper) processLookahead(in float64) float64 {
	// Read delayed sample (from laLen samples ago)
	out := c.laDelay[c.laWrite]
	// Write new input into delay line
	c.laDelay[c.laWrite] = in
	c.laWrite++
	if c.laWrite >= c.laLen {
		c.laWrite = 0
	}

	// Envelope follower on the raw input — sees peaks before they exit
	absIn := math.Abs(in)
	if absIn > c.laEnv {
		// Attack: fast rise toward the peak
		c.laEnv += c.laAttack * (absIn - c.laEnv)
	} else {
		// Release: slow decay back down
		c.laEnv += c.laRelease * (absIn - c.laEnv)
	}

	// Target gain to keep output at ceiling
	targetGain := 1.0
	if c.laEnv > c.ceiling {
		targetGain = c.ceiling / c.laEnv
	}

	// Smooth gain transitions (same attack/release as envelope)
	if targetGain < c.laGain {
		c.laGain += c.laAttack * (targetGain - c.laGain)
	} else {
		c.laGain += c.laRelease * (targetGain - c.laGain)
	}

	return out * c.laGain
}

// clip applies either soft-knee or hard clipping depending on config.
func (c *CompositeClipper) clip(x float64) float64 {
	if c.softKnee <= 0 {
		return HardClip(x, c.ceiling)
	}
	return SoftClip(x, c.ceiling, c.softKnee)
}

// SoftClip applies tanh-based soft clipping with a configurable knee.
//
// Below (ceiling - knee): linear, no distortion.
// Above (ceiling - knee): tanh compression, asymptotically approaching ceiling.
// The transition is C1-continuous (slope = 1.0 at the knee boundary).
//
// This generates significantly fewer harmonics than hard clipping, which
// means the notch filters have less work to do and produce less overshoot.
func SoftClip(x, ceiling, knee float64) float64 {
	if knee <= 0 {
		return HardClip(x, ceiling)
	}

	threshold := ceiling - knee
	if threshold < 0 {
		threshold = 0
	}

	ax := math.Abs(x)
	if ax <= threshold {
		return x // linear region — no distortion
	}

	s := 1.0
	if x < 0 {
		s = -1.0
	}

	// tanh compression: excess above threshold is compressed toward knee.
	// At excess=0: output = threshold, slope = 1.0 (C1 continuous).
	// At excess→∞: output → threshold + knee = ceiling.
	excess := ax - threshold
	compressed := threshold + knee*math.Tanh(excess/knee)
	return s * compressed
}

// Reset clears all filter and look-ahead state, as if freshly constructed.
func (c *CompositeClipper) Reset() {
	for i := range c.notch19 {
		c.notch19[i].Reset()
		c.notch57[i].Reset()
	}
	if c.laEnabled {
		for i := range c.laDelay {
			c.laDelay[i] = 0
		}
		c.laWrite = 0
		c.laEnv = 0
		c.laGain = 1.0
	}
}

// Stats returns diagnostic values for monitoring/logging.
func (c *CompositeClipper) Stats() CompositeClipperStats {
	return CompositeClipperStats{
		LookaheadGain: c.laGain,
		Envelope:      c.laEnv,
	}
}

// CompositeClipperStats exposes internal state for diagnostics.
type CompositeClipperStats struct {
	LookaheadGain float64 `json:"lookaheadGain"`
	Envelope      float64 `json:"envelope"`
}