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

package dsp

import "math"

// HilbertFilter implements a FIR-based Hilbert transform with matched delay.
// It produces the analytic signal's imaginary part (90° phase shift of all
// frequency components) while also providing the appropriately delayed
// version of the original signal.
//
// Used for SSB-SC stereo encoding (Foti/Tarsio): the L-R difference signal
// is split into in-phase (delayed) and quadrature (Hilbert-transformed)
// components for single-sideband modulation.
type HilbertFilter struct {
	coeffs []float64 // FIR coefficients (length N, odd)
	delay  []float64 // circular buffer
	pos    int       // write position in buffer
	n      int       // filter length
	mid    int       // group delay = (N-1)/2
}

// NewHilbertFilter creates a Hilbert transform FIR filter.
// taps must be odd (typical: 127 for 228 kHz). Higher = better low-freq response.
func NewHilbertFilter(taps int) *HilbertFilter {
	if taps%2 == 0 {
		taps++ // must be odd
	}
	mid := (taps - 1) / 2
	coeffs := make([]float64, taps)

	// Hilbert FIR: h[n] = 2/(π·(n-mid)) for odd (n-mid), 0 for even, Hann windowed.
	for i := 0; i < taps; i++ {
		k := i - mid
		if k == 0 {
			coeffs[i] = 0 // center tap is always zero
		} else if k%2 != 0 {
			// Odd offset: ideal Hilbert coefficient × Hann window
			window := 0.5 * (1.0 - math.Cos(2*math.Pi*float64(i)/float64(taps)))
			coeffs[i] = (2.0 / (math.Pi * float64(k))) * window
		}
		// Even offset: coefficient is zero (already initialized)
	}

	return &HilbertFilter{
		coeffs: coeffs,
		delay:  make([]float64, taps),
		pos:    0,
		n:      taps,
		mid:    mid,
	}
}

// Process takes one input sample and returns (delayed, hilbert).
// - delayed: input delayed by (N-1)/2 samples (group delay matched)
// - hilbert: 90° phase-shifted version of the input
//
// Both outputs are time-aligned: delayed[n] corresponds to hilbert[n].
func (h *HilbertFilter) Process(in float64) (delayed, hilbert float64) {
	// Write new sample
	h.delay[h.pos] = in

	// Delayed output: sample from (N-1)/2 ago
	delayIdx := (h.pos - h.mid + h.n) % h.n
	delayed = h.delay[delayIdx]

	// Hilbert output: FIR convolution
	var sum float64
	for i := 0; i < h.n; i++ {
		idx := (h.pos - i + h.n) % h.n
		sum += h.delay[idx] * h.coeffs[i]
	}
	hilbert = sum

	h.pos = (h.pos + 1) % h.n
	return
}

// Reset clears the filter state.
func (h *HilbertFilter) Reset() {
	for i := range h.delay {
		h.delay[i] = 0
	}
	h.pos = 0
}