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

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 {
	b0, b1, b2 float64
	a1, a2     float64
	z1, z2     float64 // state (direct form II transposed)
}

// 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}
	}

	omega := 2 * math.Pi * cutoffHz / sampleRate
	cosW := math.Cos(omega)
	sinW := math.Sin(omega)
	alpha := sinW / (2 * math.Sqrt2) // Q = 1/√2 for Butterworth

	a0 := 1 + alpha
	return &BiquadLPF{
		b0: (1 - cosW) / 2 / a0,
		b1: (1 - cosW) / a0,
		b2: (1 - cosW) / 2 / a0,
		a1: (-2 * cosW) / a0,
		a2: (1 - alpha) / a0,
	}
}

// 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
}

// Reset clears the filter state.
func (f *BiquadLPF) Reset() {
	f.z1 = 0
	f.z2 = 0
}