Bladeren bron

feat: add iterative composite clipper for MPX processing

main
Jan 1 maand geleden
bovenliggende
commit
1ae8b64ac8
4 gewijzigde bestanden met toevoegingen van 570 en 4 verwijderingen
  1. +30
    -1
      internal/config/config.go
  2. +311
    -0
      internal/dsp/compositeclipper.go
  3. +207
    -0
      internal/dsp/compositeclipper_test.go
  4. +22
    -3
      internal/offline/generator.go

+ 30
- 1
internal/config/config.go Bestand weergeven

@@ -49,7 +49,19 @@ type FMConfig struct {
FMModulationEnabled bool `json:"fmModulationEnabled"`
MpxGain float64 `json:"mpxGain"` // hardware calibration: scales entire composite output (default 1.0)
BS412Enabled bool `json:"bs412Enabled"` // ITU-R BS.412 MPX power limiter (EU requirement)
BS412ThresholdDBr float64 `json:"bs412ThresholdDBr"` // power limit in dBr (0 = standard, +3 = relaxed)
BS412ThresholdDBr float64 `json:"bs412ThresholdDBr"` // power limit in dBr (0 = standard, +3 = relaxed)
CompositeClipper CompositeClipperConfig `json:"compositeClipper"` // ITU-R SM.1268 iterative composite clipper
}

// CompositeClipperConfig controls the broadcast-grade composite MPX clipper.
// When enabled, replaces the simple clip→notch chain with an iterative
// clip-filter-clip processor, optionally with soft-knee clipping and
// look-ahead peak limiting.
type CompositeClipperConfig struct {
Enabled bool `json:"enabled"` // false = legacy single-pass clip
Iterations int `json:"iterations"` // clip-filter-clip passes (1-5, default 3)
SoftKnee float64 `json:"softKnee"` // soft-clip transition zone (0=hard, 0.15=moderate, 0.3=gentle)
LookaheadMs float64 `json:"lookaheadMs"` // look-ahead delay (0=off, 1.0=typical)
}

type BackendConfig struct {
@@ -157,6 +169,12 @@ func Default() Config {
LimiterCeiling: 1.0,
FMModulationEnabled: true,
MpxGain: 1.0,
CompositeClipper: CompositeClipperConfig{
Enabled: false,
Iterations: 3,
SoftKnee: 0.15,
LookaheadMs: 1.0,
},
},
Backend: BackendConfig{Kind: "file", OutputPath: "build/out/composite.f32"},
Control: ControlConfig{ListenAddress: "127.0.0.1:8088"},
@@ -319,6 +337,17 @@ func (c Config) Validate() error {
if c.FM.MpxGain < 0.1 || c.FM.MpxGain > 5 {
return fmt.Errorf("fm.mpxGain out of range (0.1..5)")
}
if c.FM.CompositeClipper.Enabled {
if c.FM.CompositeClipper.Iterations < 1 || c.FM.CompositeClipper.Iterations > 5 {
return fmt.Errorf("fm.compositeClipper.iterations out of range (1-5)")
}
if c.FM.CompositeClipper.SoftKnee < 0 || c.FM.CompositeClipper.SoftKnee > 0.5 {
return fmt.Errorf("fm.compositeClipper.softKnee out of range (0-0.5)")
}
if c.FM.CompositeClipper.LookaheadMs < 0 || c.FM.CompositeClipper.LookaheadMs > 5 {
return fmt.Errorf("fm.compositeClipper.lookaheadMs out of range (0-5)")
}
}
if c.Backend.Kind == "" {
return fmt.Errorf("backend.kind is required")
}


+ 311
- 0
internal/dsp/compositeclipper.go Bestand weergeven

@@ -0,0 +1,311 @@
// 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"`
}

+ 207
- 0
internal/dsp/compositeclipper_test.go Bestand weergeven

@@ -0,0 +1,207 @@
package dsp

import (
"math"
"testing"
)

func TestCompositeClipperPeakLimit(t *testing.T) {
// A signal at 2× ceiling must never exceed ceiling at output.
c := NewCompositeClipper(CompositeClipperConfig{
Ceiling: 1.0,
Iterations: 3,
SoftKnee: 0,
SampleRate: 228000,
})

rate := 228000.0
n := int(rate * 0.05) // 50ms
maxOut := 0.0
for i := 0; i < n; i++ {
// 1kHz tone at amplitude 2.0 (200% modulation)
in := 2.0 * math.Sin(2*math.Pi*1000*float64(i)/rate)
out := c.Process(in)
if math.Abs(out) > maxOut {
maxOut = math.Abs(out)
}
}

if maxOut > 1.001 {
t.Errorf("peak %.6f exceeds ceiling 1.0", maxOut)
}
t.Logf("hard clip: peak=%.6f (ceiling=1.0)", maxOut)
}

func TestCompositeClipperSoftKneePeakLimit(t *testing.T) {
c := NewCompositeClipper(CompositeClipperConfig{
Ceiling: 1.0,
Iterations: 3,
SoftKnee: 0.2,
SampleRate: 228000,
})

rate := 228000.0
n := int(rate * 0.05)
maxOut := 0.0
for i := 0; i < n; i++ {
in := 2.0 * math.Sin(2*math.Pi*1000*float64(i)/rate)
out := c.Process(in)
if math.Abs(out) > maxOut {
maxOut = math.Abs(out)
}
}

if maxOut > 1.001 {
t.Errorf("soft clip peak %.6f exceeds ceiling 1.0", maxOut)
}
t.Logf("soft clip (knee=0.2): peak=%.6f", maxOut)
}

func TestCompositeClipperLookaheadPeakLimit(t *testing.T) {
c := NewCompositeClipper(CompositeClipperConfig{
Ceiling: 1.0,
Iterations: 3,
SoftKnee: 0.15,
LookaheadMs: 1.0,
SampleRate: 228000,
})

rate := 228000.0
n := int(rate * 0.10) // 100ms to let look-ahead settle
maxOut := 0.0
for i := 0; i < n; i++ {
in := 2.0 * math.Sin(2*math.Pi*1000*float64(i)/rate)
out := c.Process(in)
if math.Abs(out) > maxOut {
maxOut = math.Abs(out)
}
}

if maxOut > 1.001 {
t.Errorf("lookahead peak %.6f exceeds ceiling 1.0", maxOut)
}
t.Logf("lookahead (1ms, soft knee=0.15, 3 iter): peak=%.6f", maxOut)
}

func TestCompositeClipperGuardBands(t *testing.T) {
// Verify that after clipping, 19kHz and 57kHz energy is suppressed.
// Feed a 1kHz sine at high level → clips → generates harmonics.
// The clipper's notch filters should remove 19kHz and 57kHz.
rate := 228000.0
n := int(rate * 0.1) // 100ms

for _, tc := range []struct {
name string
iterations int
softKnee float64
lookahead float64
}{
{"hard-1iter", 1, 0, 0},
{"hard-3iter", 3, 0, 0},
{"soft-3iter", 3, 0.15, 0},
{"lookahead-3iter", 3, 0.15, 1.0},
} {
t.Run(tc.name, func(t *testing.T) {
c := NewCompositeClipper(CompositeClipperConfig{
Ceiling: 1.0,
Iterations: tc.iterations,
SoftKnee: tc.softKnee,
LookaheadMs: tc.lookahead,
SampleRate: rate,
})

out := make([]float64, n)
for i := 0; i < n; i++ {
in := 2.0 * math.Sin(2*math.Pi*1000*float64(i)/rate)
out[i] = c.Process(in)
}

// Skip first 10ms for filter settling
skip := int(rate * 0.01)
analysis := out[skip:]

e19 := GoertzelEnergy(analysis, rate, 19000)
e57 := GoertzelEnergy(analysis, rate, 57000)
e1k := GoertzelEnergy(analysis, rate, 1000)

r19 := -100.0
r57 := -100.0
if e1k > 0 {
if e19 > 0 {
r19 = 10 * math.Log10(e19/e1k)
}
if e57 > 0 {
r57 = 10 * math.Log10(e57/e1k)
}
}

t.Logf("19kHz: %.1f dB below 1kHz, 57kHz: %.1f dB below 1kHz", -r19, -r57)

// With 3 iterations, 19kHz should be suppressed by at least 40dB
if tc.iterations >= 3 && r19 > -40 {
t.Errorf("19kHz not sufficiently suppressed: %.1f dB (want < -40 dB)", r19)
}
})
}
}

func TestSoftClipContinuity(t *testing.T) {
// Verify C1 continuity at the knee boundary
ceiling := 1.0
knee := 0.2
threshold := ceiling - knee // 0.8

// Slope just below threshold (linear region)
dx := 1e-8
y0 := SoftClip(threshold-dx, ceiling, knee)
y1 := SoftClip(threshold+dx, ceiling, knee)
slope := (y1 - y0) / (2 * dx)

if math.Abs(slope-1.0) > 0.01 {
t.Errorf("slope at knee boundary = %.4f, want ~1.0", slope)
}

// Verify output never exceeds ceiling
for x := 0.0; x <= 5.0; x += 0.001 {
y := SoftClip(x, ceiling, knee)
if y > ceiling+1e-9 {
t.Fatalf("SoftClip(%.3f) = %.6f > ceiling %.6f", x, y, ceiling)
}
}

// Verify asymptotic approach to ceiling
y5 := SoftClip(5.0, ceiling, knee)
if y5 < 0.99 {
t.Errorf("SoftClip(5.0) = %.6f, expected near ceiling (%.6f)", y5, ceiling)
}
}

func TestCompositeClipperReset(t *testing.T) {
c := NewCompositeClipper(CompositeClipperConfig{
Ceiling: 1.0,
Iterations: 2,
SoftKnee: 0.15,
LookaheadMs: 1.0,
SampleRate: 228000,
})

// Process some samples
for i := 0; i < 1000; i++ {
c.Process(0.5)
}

stats := c.Stats()
if stats.Envelope == 0 {
t.Error("envelope should be non-zero after processing")
}

c.Reset()

stats = c.Stats()
if stats.Envelope != 0 {
t.Errorf("envelope should be 0 after reset, got %f", stats.Envelope)
}
if stats.LookaheadGain != 1.0 {
t.Errorf("gain should be 1.0 after reset, got %f", stats.LookaheadGain)
}
}

+ 22
- 3
internal/offline/generator.go Bestand weergeven

@@ -41,6 +41,8 @@ type LiveParams struct {
ToneRightHz float64
ToneAmplitude float64
AudioGain float64
// Composite clipper: live-toggleable without DSP chain reinit.
CompositeClipperEnabled bool
}

// PreEmphasizedSource wraps an audio source and applies pre-emphasis.
@@ -110,6 +112,7 @@ type Generator struct {
mpxNotch19 *dsp.FilterChain // composite clipper protection
mpxNotch57 *dsp.FilterChain
bs412 *dsp.BS412Limiter // ITU-R BS.412 MPX power limiter (optional)
compositeClip *dsp.CompositeClipper // ITU-R SM.1268 iterative composite clipper (optional)

// Pre-allocated frame buffer — reused every GenerateFrame call.
frameBuf *output.CompositeFrame
@@ -223,6 +226,15 @@ func (g *Generator) init() {
g.cleanupLPF_R = dsp.NewAudioLPF(g.sampleRate)
// Composite clipper protection: double-notch at 19kHz + 57kHz
g.mpxNotch19, g.mpxNotch57 = dsp.NewCompositeProtection(g.sampleRate)
// ITU-R SM.1268 iterative composite clipper (optional, replaces simple clip+notch)
// Always created so it can be live-toggled via CompositeClipperEnabled.
g.compositeClip = dsp.NewCompositeClipper(dsp.CompositeClipperConfig{
Ceiling: ceiling,
Iterations: g.cfg.FM.CompositeClipper.Iterations,
SoftKnee: g.cfg.FM.CompositeClipper.SoftKnee,
LookaheadMs: g.cfg.FM.CompositeClipper.LookaheadMs,
SampleRate: g.sampleRate,
})
// BS.412 MPX power limiter (EU/CH requirement for licensed FM)
if g.cfg.FM.BS412Enabled {
// chunkSec is not known at init time (Engine.chunkDuration may differ).
@@ -260,6 +272,7 @@ func (g *Generator) init() {
ToneRightHz: g.cfg.Audio.ToneRightHz,
ToneAmplitude: g.cfg.Audio.ToneAmplitude,
AudioGain: g.cfg.Audio.Gain,
CompositeClipperEnabled: g.cfg.FM.CompositeClipper.Enabled,
})

if g.licenseState != nil {
@@ -418,9 +431,15 @@ func (g *Generator) GenerateFrame(duration time.Duration) *output.CompositeFrame
if lp.StereoEnabled {
audioMPX += float64(comps.Stereo)
}
audioMPX = dsp.HardClip(audioMPX, ceiling)
audioMPX = g.mpxNotch19.Process(audioMPX)
audioMPX = g.mpxNotch57.Process(audioMPX)
if lp.CompositeClipperEnabled && g.compositeClip != nil {
// ITU-R SM.1268 iterative clipper: look-ahead + N×(clip→notch→notch) + final clip
audioMPX = g.compositeClip.Process(audioMPX)
} else {
// Legacy single-pass: one clip, then notch, no final safety clip
audioMPX = dsp.HardClip(audioMPX, ceiling)
audioMPX = g.mpxNotch19.Process(audioMPX)
audioMPX = g.mpxNotch57.Process(audioMPX)
}

// BS.412: apply gain and measure power
if bs412Gain < 1.0 {


Laden…
Annuleren
Opslaan