refactor: harden broadcast DSP chain around pilot and RDS protection

Upgrade the audio path to a broadcast-style chain with a 4th-order 15 kHz low-pass, 19 kHz pilot notch, audio-only composite clipping, and post-clip protection notches at 19/57 kHz before adding clean pilot and RDS carriers. This reorders the processing so pilot and RDS stay fixed and unclipped while audio dynamics are constrained within the available modulation budget.

internal/dsp/biquad.go

@@ -2,35 +2,63 @@ package dsp
 
 import "math"
 
-// BiquadLPF is a second-order Butterworth lowpass filter (biquad, direct form II transposed).
-// Used after the audio limiter to remove intermodulation products and harmonics
-// that could fall into the 19kHz pilot, 38kHz stereo sub, or 57kHz RDS bands.
-//
-// At 228kHz with fc=15kHz:
-//   15kHz: -3 dB (corner)
-//   19kHz: -5 dB
-//   38kHz: -18 dB
-//   57kHz: -27 dB  ← protects RDS band
-type BiquadLPF struct {
+// Biquad is a generic second-order IIR filter (direct form II transposed).
+type Biquad struct {
 	b0, b1, b2 float64
 	a1, a2     float64
-	z1, z2     float64 // state (direct form II transposed)
+	z1, z2     float64
 }
 
-// NewBiquadLPF creates a 2nd-order Butterworth lowpass at the given cutoff.
-func NewBiquadLPF(cutoffHz, sampleRate float64) *BiquadLPF {
-	if cutoffHz <= 0 || sampleRate <= 0 || cutoffHz >= sampleRate/2 {
-		// Passthrough: return unity filter
-		return &BiquadLPF{b0: 1}
+// Process filters one sample.
+func (f *Biquad) Process(in float64) float64 {
+	out := f.b0*in + f.z1
+	f.z1 = f.b1*in - f.a1*out + f.z2
+	f.z2 = f.b2*in - f.a2*out
+	return out
+}
+
+// Reset clears state.
+func (f *Biquad) Reset() { f.z1 = 0; f.z2 = 0 }
+
+// FilterChain cascades multiple biquad sections in series.
+// Used for higher-order filters (e.g. 4th-order = 2 biquads).
+type FilterChain struct {
+	Stages []Biquad
+}
+
+// Process runs input through all stages in series.
+func (c *FilterChain) Process(in float64) float64 {
+	x := in
+	for i := range c.Stages {
+		x = c.Stages[i].Process(x)
 	}
+	return x
+}
+
+// Reset clears all filter state.
+func (c *FilterChain) Reset() {
+	for i := range c.Stages {
+		c.Stages[i].Reset()
+	}
+}
 
+// --- Factory functions ---
+
+// NewBiquadLPF creates a 2nd-order Butterworth lowpass (Q = 1/√2).
+func NewBiquadLPF(cutoffHz, sampleRate float64) *Biquad {
+	return newBiquadLPFWithQ(cutoffHz, sampleRate, math.Sqrt2/2)
+}
+
+func newBiquadLPFWithQ(cutoffHz, sampleRate, q float64) *Biquad {
+	if cutoffHz <= 0 || sampleRate <= 0 || cutoffHz >= sampleRate/2 {
+		return &Biquad{b0: 1} // passthrough
+	}
 	omega := 2 * math.Pi * cutoffHz / sampleRate
 	cosW := math.Cos(omega)
 	sinW := math.Sin(omega)
-	alpha := sinW / (2 * math.Sqrt2) // Q = 1/√2 for Butterworth
-
+	alpha := sinW / (2 * q)
 	a0 := 1 + alpha
-	return &BiquadLPF{
+	return &Biquad{
 		b0: (1 - cosW) / 2 / a0,
 		b1: (1 - cosW) / a0,
 		b2: (1 - cosW) / 2 / a0,
@@ -39,16 +67,64 @@ func NewBiquadLPF(cutoffHz, sampleRate float64) *BiquadLPF {
 	}
 }
 
-// Process filters a single sample.
-func (f *BiquadLPF) Process(in float64) float64 {
-	out := f.b0*in + f.z1
-	f.z1 = f.b1*in - f.a1*out + f.z2
-	f.z2 = f.b2*in - f.a2*out
-	return out
+// NewLPF4 creates a 4th-order Butterworth lowpass (two cascaded biquads).
+// Provides -24dB/octave rolloff. At 228kHz with fc=15kHz:
+//
+//	15kHz: -6dB, 19kHz: -14dB, 38kHz: -36dB, 57kHz: -54dB
+func NewLPF4(cutoffHz, sampleRate float64) *FilterChain {
+	// 4th-order Butterworth: cascade two 2nd-order sections with Q values
+	// derived from the pole angles: π/8 and 3π/8
+	q1 := 1.0 / (2 * math.Cos(math.Pi/8))   // ≈ 0.5412
+	q2 := 1.0 / (2 * math.Cos(3*math.Pi/8))  // ≈ 1.3066
+	return &FilterChain{
+		Stages: []Biquad{
+			*newBiquadLPFWithQ(cutoffHz, sampleRate, q1),
+			*newBiquadLPFWithQ(cutoffHz, sampleRate, q2),
+		},
+	}
+}
+
+// NewNotch creates a 2nd-order IIR notch (bandstop) filter.
+// Q controls width: higher Q = narrower notch.
+// Typical: Q=5 → ~4kHz wide at -3dB, Q=10 → ~2kHz wide.
+func NewNotch(centerHz, sampleRate, q float64) *Biquad {
+	if centerHz <= 0 || sampleRate <= 0 || centerHz >= sampleRate/2 {
+		return &Biquad{b0: 1}
+	}
+	omega := 2 * math.Pi * centerHz / sampleRate
+	cosW := math.Cos(omega)
+	alpha := math.Sin(omega) / (2 * q)
+	a0 := 1 + alpha
+	return &Biquad{
+		b0: 1 / a0,
+		b1: -2 * cosW / a0,
+		b2: 1 / a0,
+		a1: -2 * cosW / a0,
+		a2: (1 - alpha) / a0,
+	}
+}
+
+// --- Broadcast-specific filter factories ---
+
+// NewAudioLPF creates the broadcast-standard 15kHz audio lowpass.
+// 4th-order Butterworth ensures the guard band (15–23kHz) is clean.
+func NewAudioLPF(sampleRate float64) *FilterChain {
+	return NewLPF4(15000, sampleRate)
+}
+
+// NewPilotNotch creates a narrow notch at 19kHz to kill residual audio
+// energy at the pilot frequency. Applied BEFORE stereo encoding.
+// Q=5: -3dB width ~4kHz, >40dB rejection at 19kHz, <0.5dB loss at 15kHz.
+func NewPilotNotch(sampleRate float64) *Biquad {
+	return NewNotch(19000, sampleRate, 5)
 }
 
-// Reset clears the filter state.
-func (f *BiquadLPF) Reset() {
-	f.z1 = 0
-	f.z2 = 0
+// NewCompositeProtection creates notch filters for the composite clipper.
+// Applied to clipped audio composite (mono+stereo_sub) to remove clip
+// harmonics from the pilot (19kHz) and RDS (57kHz) bands before adding
+// the actual pilot and RDS signals at clean, fixed levels.
+func NewCompositeProtection(sampleRate float64) (notch19, notch57 *Biquad) {
+	notch19 = NewNotch(19000, sampleRate, 3) // wider Q for broadband clip artifacts
+	notch57 = NewNotch(57000, sampleRate, 3)
+	return
 }

internal/offline/generator.go

@@ -77,13 +77,24 @@ type Generator struct {
 	stereoEncoder stereo.StereoEncoder
 	rdsEnc        *rds.Encoder
 	combiner      mpx.DefaultCombiner
-	limiter       *dsp.StereoLimiter // stereo-linked, operates on L/R BEFORE stereo encoding
-	lpfL, lpfR    *dsp.BiquadLPF    // 15kHz lowpass after limiter, protects RDS band
+	limiter       *dsp.StereoLimiter  // stereo-linked, operates on L/R BEFORE stereo encoding
 	fmMod         *dsp.FMModulator
 	sampleRate    float64
 	initialized   bool
 	frameSeq      uint64
 
+	// Broadcast-standard audio filter chain (per channel, L and R):
+	//   Pre-emphasis → 15kHz LPF (4th-order) → 19kHz Notch → Drive → Limiter
+	audioLPF_L    *dsp.FilterChain  // 4th-order Butterworth 15kHz
+	audioLPF_R    *dsp.FilterChain
+	pilotNotchL   *dsp.Biquad       // 19kHz notch (guard band protection)
+	pilotNotchR   *dsp.Biquad
+
+	// Composite clipper protection (post-clip notch filters):
+	//   Audio composite → clip → notch 19kHz → notch 57kHz → + pilot → + RDS
+	mpxNotch19    *dsp.Biquad       // removes clip harmonics at pilot freq
+	mpxNotch57    *dsp.Biquad       // removes clip harmonics at RDS freq
+
 	// Pre-allocated frame buffer — reused every GenerateFrame call.
 	frameBuf *output.CompositeFrame
 	bufCap   int
@@ -160,10 +171,16 @@ func (g *Generator) init() {
 	if g.cfg.FM.LimiterEnabled {
 		g.limiter = dsp.NewStereoLimiter(audioCeiling, 0.5, 100, g.sampleRate)
 	}
-	// 15kHz lowpass after limiter — removes limiter gain-step intermodulation
-	// products that would otherwise fall into pilot/stereo/RDS bands.
-	g.lpfL = dsp.NewBiquadLPF(15000, g.sampleRate)
-	g.lpfR = dsp.NewBiquadLPF(15000, g.sampleRate)
+
+	// Broadcast-standard filter chain:
+	// 1) 15kHz 4th-order Butterworth LPF — steep guard band, -14dB@19kHz, -54dB@57kHz
+	g.audioLPF_L = dsp.NewAudioLPF(g.sampleRate)
+	g.audioLPF_R = dsp.NewAudioLPF(g.sampleRate)
+	// 2) 19kHz notch — kills residual audio at pilot freq, >40dB rejection
+	g.pilotNotchL = dsp.NewPilotNotch(g.sampleRate)
+	g.pilotNotchR = dsp.NewPilotNotch(g.sampleRate)
+	// 3) Composite clipper protection notches at 19kHz + 57kHz
+	g.mpxNotch19, g.mpxNotch57 = dsp.NewCompositeProtection(g.sampleRate)
 	if g.cfg.FM.FMModulationEnabled {
 		g.fmMod = dsp.NewFMModulator(g.sampleRate)
 		if g.cfg.FM.MaxDeviationHz > 0 { g.fmMod.MaxDeviation = g.cfg.FM.MaxDeviationHz }
@@ -223,52 +240,71 @@ func (g *Generator) GenerateFrame(duration time.Duration) *output.CompositeFrame
 		lp = &LiveParams{OutputDrive: 1.0, LimiterCeiling: 1.0}
 	}
 
-	// Signal path (matches professional broadcast processors):
-	//   Audio L/R → × Drive → Stereo-linked limiter → Stereo encoder
-	//   → Mono + Stereo sub (from limited audio, natural levels)
-	//   → + Pilot (fixed) → + RDS (fixed) → FM modulator
+	// Broadcast-standard FM MPX signal chain:
 	//
-	// The limiter never sees the 38kHz subcarrier, so it can't pump
-	// the stereo difference signal. Pilot and RDS are post-encoder
-	// at fixed amplitudes, unaffected by audio dynamics.
+	//   Audio L/R
+	//     → PreEmphasis (50µs EU / 75µs US)
+	//     → 15kHz LPF (4th-order Butterworth, -14dB@19kHz, -54dB@57kHz)
+	//     → 19kHz Notch (>40dB rejection, guard band protection)
+	//     → × OutputDrive
+	//     → StereoLimiter (instant attack, smooth release)
+	//     → Stereo Encode → Mono (L+R)/2 + Stereo Sub (L-R)/2 × 38kHz
+	//   Audio MPX composite
+	//     → HardClip at audioCeiling (catches limiter overshoots)
+	//     → 19kHz Notch (removes clip harmonics at pilot freq)
+	//     → 57kHz Notch (removes clip harmonics at RDS freq)
+	//   + Pilot 19kHz (fixed amplitude, post-clip)
+	//   + RDS 57kHz (fixed amplitude, post-clip)
+	//   → FM Modulator
 	//
-	// Audio ceiling is auto-reduced to leave headroom for pilot + RDS,
-	// so total composite stays within ±ceiling (= ±75kHz deviation).
+	// Key: Pilot and RDS are NEVER clipped or filtered. They're added
+	// after all audio processing at constant amplitude.
 	ceiling := lp.LimiterCeiling
 	if ceiling <= 0 { ceiling = 1.0 }
 	pilotAmp := lp.PilotLevel * lp.OutputDrive
 	rdsAmp := lp.RDSInjection * lp.OutputDrive
 	audioCeiling := ceiling - pilotAmp - rdsAmp
-	if audioCeiling < 0.3 { audioCeiling = 0.3 } // safety floor
+	if audioCeiling < 0.3 { audioCeiling = 0.3 }
 
 	for i := 0; i < samples; i++ {
 		in := g.source.NextFrame()
 
-		// --- Stage 1: Band-limit pre-emphasized audio ---
-		// The 15kHz LPF goes BEFORE drive+limiter. Pre-emphasis boosts
-		// HF by up to +13.5dB. Without the LPF, the limiter would waste
-		// gain reduction on HF peaks that get filtered later, causing
-		// wild modulation swings (30-163%). With LPF first, the limiter
-		// sees the final audio bandwidth and sets gain correctly.
-		l := g.lpfL.Process(float64(in.L))
-		r := g.lpfR.Process(float64(in.R))
-
-		// --- Stage 2: Scale and limit ---
+		// --- Stage 1: Audio filtering (per-channel) ---
+		// 15kHz LPF removes out-of-band pre-emphasis energy.
+		// 19kHz notch kills residual energy at pilot frequency.
+		// Both run BEFORE drive+limiter so the limiter sees the
+		// actual audio bandwidth, not wasted HF energy.
+		l := g.audioLPF_L.Process(float64(in.L))
+		l = g.pilotNotchL.Process(l)
+		r := g.audioLPF_R.Process(float64(in.R))
+		r = g.pilotNotchR.Process(r)
+
+		// --- Stage 2: Drive + Limit ---
 		l *= lp.OutputDrive
 		r *= lp.OutputDrive
-
 		if lp.LimiterEnabled && g.limiter != nil {
 			l, r = g.limiter.Process(l, r)
 		}
 
-		// --- Stage 3: Stereo encode the limited, filtered audio ---
+		// --- Stage 3: Stereo encode ---
 		limited := audio.NewFrame(audio.Sample(l), audio.Sample(r))
 		comps := g.stereoEncoder.Encode(limited)
 
-		// --- Stage 3: Combine at fixed levels ---
-		composite := float64(comps.Mono)
+		// --- Stage 4: Audio composite clip + protection ---
+		// Clip the audio-only composite (mono + stereo sub) to budget.
+		// Then notch-filter the clip harmonics out of the pilot (19kHz)
+		// and RDS (57kHz) bands before adding the real pilot and RDS.
+		audioMPX := float64(comps.Mono)
+		if lp.StereoEnabled {
+			audioMPX += float64(comps.Stereo)
+		}
+		audioMPX = dsp.HardClip(audioMPX, audioCeiling)
+		audioMPX = g.mpxNotch19.Process(audioMPX)
+		audioMPX = g.mpxNotch57.Process(audioMPX)
+
+		// --- Stage 5: Add protected components at fixed levels ---
+		composite := audioMPX
 		if lp.StereoEnabled {
-			composite += float64(comps.Stereo)
 			composite += pilotAmp * comps.Pilot
 		}
 		if g.rdsEnc != nil && lp.RDSEnabled {
@@ -277,13 +313,6 @@ func (g *Generator) GenerateFrame(duration time.Duration) *output.CompositeFrame
 			composite += rdsAmp * rdsValue
 		}
 
-		// Final composite safety clip — only fires on brief limiter
-		// overshoots during fast transients. Clips the entire composite,
-		// not individual audio bands, so harmonics don't target RDS.
-		if lp.LimiterEnabled {
-			composite = dsp.HardClip(composite, ceiling)
-		}
-
 		if g.fmMod != nil {
 			iq_i, iq_q := g.fmMod.Modulate(composite)
 			frame.Samples[i] = output.IQSample{I: float32(iq_i), Q: float32(iq_q)}