From c22fb7bf2a3162b3bdb5da2239a8f2ed5f8d3f8c Mon Sep 17 00:00:00 2001 From: Jan Svabenik Date: Fri, 3 Apr 2026 09:13:15 +0200 Subject: [PATCH] feat: add pre-emphasis, MPX limiter and FM IQ modulation DSP blocks --- go.mod | 2 +- internal/dsp/fmmod.go | 46 ++++++++++++++++ internal/dsp/fmmod_test.go | 41 ++++++++++++++ internal/dsp/limiter.go | 70 ++++++++++++++++++++++++ internal/dsp/limiter_test.go | 46 ++++++++++++++++ internal/dsp/preemphasis.go | 94 ++++++++++++++++++++++++++++++++ internal/dsp/preemphasis_test.go | 72 ++++++++++++++++++++++++ 7 files changed, 370 insertions(+), 1 deletion(-) create mode 100644 internal/dsp/fmmod.go create mode 100644 internal/dsp/fmmod_test.go create mode 100644 internal/dsp/limiter.go create mode 100644 internal/dsp/limiter_test.go create mode 100644 internal/dsp/preemphasis.go create mode 100644 internal/dsp/preemphasis_test.go diff --git a/go.mod b/go.mod index 568d50e..57c4f25 100644 --- a/go.mod +++ b/go.mod @@ -1,6 +1,6 @@ module github.com/jan/fm-rds-tx -go 1.24.0 +go 1.22 require github.com/jan/fm-rds-tx/internal v0.0.0 diff --git a/internal/dsp/fmmod.go b/internal/dsp/fmmod.go new file mode 100644 index 0000000..9036f2a --- /dev/null +++ b/internal/dsp/fmmod.go @@ -0,0 +1,46 @@ +package dsp + +import "math" + +// FMModulator converts a composite baseband signal into FM-modulated IQ samples. +// The modulated output represents the instantaneous frequency deviation around +// a zero-IF carrier (baseband IQ). For a real SDR transmitter the SDR hardware +// then upconverts this to the desired center frequency. +type FMModulator struct { + // MaxDeviation is the peak frequency deviation in Hz (±75 kHz for FM broadcast). + MaxDeviation float64 + SampleRate float64 + + phase float64 // accumulated carrier phase in radians +} + +// NewFMModulator creates a modulator with broadcast FM defaults. +func NewFMModulator(sampleRate float64) *FMModulator { + return &FMModulator{ + MaxDeviation: 75000, // ±75 kHz + SampleRate: sampleRate, + } +} + +// Modulate converts a single composite sample (normalized to [-1,+1] representing +// full deviation) into an IQ pair. +func (m *FMModulator) Modulate(composite float64) (i, q float64) { + // Instantaneous frequency offset = composite * maxDeviation + // Phase increment per sample = 2π * freq_offset / sampleRate + freqOffset := composite * m.MaxDeviation + m.phase += 2 * math.Pi * freqOffset / m.SampleRate + + // Keep phase bounded to avoid float64 precision loss over long runs + if m.phase > math.Pi { + m.phase -= 2 * math.Pi * math.Floor((m.phase+math.Pi)/(2*math.Pi)) + } + + i = math.Cos(m.phase) + q = math.Sin(m.phase) + return i, q +} + +// Reset clears the modulator phase. +func (m *FMModulator) Reset() { + m.phase = 0 +} diff --git a/internal/dsp/fmmod_test.go b/internal/dsp/fmmod_test.go new file mode 100644 index 0000000..06f48f8 --- /dev/null +++ b/internal/dsp/fmmod_test.go @@ -0,0 +1,41 @@ +package dsp + +import ( + "math" + "testing" +) + +func TestFMModulatorOutputMagnitude(t *testing.T) { + mod := NewFMModulator(228000) + // IQ output should always have magnitude ~1 + for n := 0; n < 1000; n++ { + composite := math.Sin(2 * math.Pi * 1000 * float64(n) / 228000) + i, q := mod.Modulate(composite) + mag := math.Sqrt(i*i + q*q) + if math.Abs(mag-1.0) > 1e-9 { + t.Fatalf("sample %d: magnitude=%.9f, expected 1.0", n, mag) + } + } +} + +func TestFMModulatorZeroInput(t *testing.T) { + mod := NewFMModulator(228000) + // Zero input = no deviation = phase stays at 0 + for n := 0; n < 100; n++ { + i, q := mod.Modulate(0) + if math.Abs(i-1.0) > 1e-9 || math.Abs(q) > 1e-9 { + t.Fatalf("sample %d: expected (1,0) got (%.6f,%.6f)", n, i, q) + } + } +} + +func TestFMModulatorReset(t *testing.T) { + mod := NewFMModulator(228000) + mod.Modulate(0.5) + mod.Modulate(0.5) + mod.Reset() + i, q := mod.Modulate(0) + if math.Abs(i-1.0) > 1e-9 || math.Abs(q) > 1e-9 { + t.Fatalf("after reset: expected (1,0) got (%.6f,%.6f)", i, q) + } +} diff --git a/internal/dsp/limiter.go b/internal/dsp/limiter.go new file mode 100644 index 0000000..7fc72b4 --- /dev/null +++ b/internal/dsp/limiter.go @@ -0,0 +1,70 @@ +package dsp + +import "math" + +// MPXLimiter is a look-ahead peak limiter for the composite MPX signal. +// It prevents overmodulation by keeping the composite output within +// the configured ceiling. The limiter applies gain reduction smoothly +// to avoid introducing audible artifacts. +type MPXLimiter struct { + ceiling float64 // maximum absolute output level (e.g. 1.0) + attackCoeff float64 // attack smoothing coefficient + releaseCoeff float64 // release smoothing coefficient + gainReduction float64 // current gain reduction in linear scale (0 = no reduction) +} + +// NewMPXLimiter creates a limiter with the given ceiling, attack and release +// times in milliseconds, and sample rate. +func NewMPXLimiter(ceiling, attackMs, releaseMs, sampleRate float64) *MPXLimiter { + if ceiling <= 0 { + ceiling = 1.0 + } + if attackMs <= 0 { + attackMs = 0.1 // fast attack for MPX + } + if releaseMs <= 0 { + releaseMs = 50 // moderate release + } + attackSamples := attackMs * sampleRate / 1000 + releaseSamples := releaseMs * sampleRate / 1000 + + return &MPXLimiter{ + ceiling: ceiling, + attackCoeff: 1.0 - math.Exp(-1.0/attackSamples), + releaseCoeff: 1.0 - math.Exp(-1.0/releaseSamples), + } +} + +// Process applies limiting to a single composite sample. +func (l *MPXLimiter) Process(in float64) float64 { + absIn := math.Abs(in) + targetReduction := 0.0 + if absIn > l.ceiling { + targetReduction = 1.0 - l.ceiling/absIn + } + + if targetReduction > l.gainReduction { + l.gainReduction += l.attackCoeff * (targetReduction - l.gainReduction) + } else { + l.gainReduction += l.releaseCoeff * (targetReduction - l.gainReduction) + } + + gain := 1.0 - l.gainReduction + return in * gain +} + +// Reset clears the limiter state. +func (l *MPXLimiter) Reset() { + l.gainReduction = 0 +} + +// HardClip provides a simple hard clipper as a safety net after the limiter. +func HardClip(sample, ceiling float64) float64 { + if sample > ceiling { + return ceiling + } + if sample < -ceiling { + return -ceiling + } + return sample +} diff --git a/internal/dsp/limiter_test.go b/internal/dsp/limiter_test.go new file mode 100644 index 0000000..240cf75 --- /dev/null +++ b/internal/dsp/limiter_test.go @@ -0,0 +1,46 @@ +package dsp + +import ( + "math" + "testing" +) + +func TestMPXLimiterPassesQuiet(t *testing.T) { + lim := NewMPXLimiter(1.0, 0.1, 50, 228000) + for i := 0; i < 100; i++ { + in := 0.3 * math.Sin(float64(i)) + out := lim.Process(in) + if math.Abs(out-in) > 1e-9 { + t.Fatalf("limiter altered quiet signal at sample %d: in=%.6f out=%.6f", i, in, out) + } + } +} + +func TestMPXLimiterClamps(t *testing.T) { + lim := NewMPXLimiter(1.0, 0.01, 50, 228000) + // Feed a signal well above ceiling + var maxOut float64 + for i := 0; i < 10000; i++ { + in := 3.0 * math.Sin(2*math.Pi*1000*float64(i)/228000) + out := lim.Process(in) + if math.Abs(out) > maxOut { + maxOut = math.Abs(out) + } + } + // After attack settles, output should approach ceiling + if maxOut > 1.5 { + t.Fatalf("limiter didn't reduce level enough: maxOut=%.4f", maxOut) + } +} + +func TestHardClip(t *testing.T) { + if HardClip(1.5, 1.0) != 1.0 { + t.Fatal("expected clip to 1.0") + } + if HardClip(-1.5, 1.0) != -1.0 { + t.Fatal("expected clip to -1.0") + } + if HardClip(0.5, 1.0) != 0.5 { + t.Fatal("expected passthrough") + } +} diff --git a/internal/dsp/preemphasis.go b/internal/dsp/preemphasis.go new file mode 100644 index 0000000..d627aad --- /dev/null +++ b/internal/dsp/preemphasis.go @@ -0,0 +1,94 @@ +package dsp + +import "math" + +// PreEmphasis implements a first-order high-shelf filter for FM broadcast +// pre-emphasis. Standard time constants are 50 µs (Europe/world) and +// 75 µs (North America/South Korea). +// +// Transfer function: H(s) = 1 + s*τ +// Bilinear transform to discrete: H(z) = (b0 + b1*z^-1) / (1 + a1*z^-1) +type PreEmphasis struct { + b0, b1, a1 float64 + x1, y1 float64 // state + enabled bool +} + +// NewPreEmphasis creates a pre-emphasis filter for the given time constant +// and sample rate. tau is in microseconds (50 or 75). +func NewPreEmphasis(tauMicroseconds, sampleRate float64) *PreEmphasis { + if tauMicroseconds <= 0 || sampleRate <= 0 { + return &PreEmphasis{enabled: false} + } + tau := tauMicroseconds * 1e-6 + + // Bilinear transform of H(s) = 1 + s*tau + // With pre-warping: let c = 2*fs + // Numerator: (1 + tau*c) + (-1 + tau*c)*z^-1 => b0 = 1+tau*c, b1 = -1+tau*c + // Denominator: (1 + tau*c) + (1 - tau*c)*z^-1 (but we normalize so a0=1) + // Wait - cleaner approach: standard first-order shelf. + + // H(s) = (1 + s*tau) mapped with bilinear: s = c*(1 - z^-1)/(1 + z^-1), c = 2*fs + // De-emphasis: H_de(z) = (1-alpha)/(1 - alpha*z^-1), alpha = exp(-1/(tau*fs)) + // Pre-emphasis: H_pre(z) = 1/H_de(z) = (1 - alpha*z^-1)/(1-alpha) + // y[n] = (x[n] - alpha*x[n-1]) / (1 - alpha) + alpha := math.Exp(-1.0 / (tau * sampleRate)) + gain := 1.0 / (1.0 - alpha) + + return &PreEmphasis{ + b0: gain, + b1: -alpha * gain, + a1: 0, // FIR, no feedback + enabled: true, + } +} + +// Process applies the pre-emphasis filter to a single sample. +func (p *PreEmphasis) Process(in float64) float64 { + if !p.enabled { + return in + } + out := p.b0*in + p.b1*p.x1 + p.x1 = in + return out +} + +// Reset clears the filter state. +func (p *PreEmphasis) Reset() { + p.x1 = 0 + p.y1 = 0 +} + +// DeEmphasis implements the complementary de-emphasis filter. +// H(z) = (1-alpha)/(1 - alpha*z^-1) +type DeEmphasis struct { + alpha float64 + gain float64 + prevOut float64 + enabled bool +} + +// NewDeEmphasis creates a de-emphasis filter. +func NewDeEmphasis(tauMicroseconds, sampleRate float64) *DeEmphasis { + if tauMicroseconds <= 0 || sampleRate <= 0 { + return &DeEmphasis{enabled: false} + } + tau := tauMicroseconds * 1e-6 + alpha := math.Exp(-1.0 / (tau * sampleRate)) + return &DeEmphasis{alpha: alpha, gain: 1.0 - alpha, enabled: true} +} + +// Process applies the de-emphasis filter. +func (d *DeEmphasis) Process(in float64) float64 { + if !d.enabled { + return in + } + out := d.gain*in + d.alpha*d.prevOut + d.prevOut = out + return out +} + +// Reset clears the filter state. +func (d *DeEmphasis) Reset() { + d.prevOut = 0 +} diff --git a/internal/dsp/preemphasis_test.go b/internal/dsp/preemphasis_test.go new file mode 100644 index 0000000..afd369f --- /dev/null +++ b/internal/dsp/preemphasis_test.go @@ -0,0 +1,72 @@ +package dsp + +import ( + "math" + "testing" +) + +func TestPreEmphasisBoostsHighFrequency(t *testing.T) { + fs := 48000.0 + pe := NewPreEmphasis(50, fs) + + // Generate a low-frequency tone and a high-frequency tone, + // measure the energy ratio after pre-emphasis. + n := 4800 // 100ms + lowFreq := 100.0 + highFreq := 10000.0 + + var lowEnergy, highEnergy float64 + for i := 0; i < n; i++ { + s := math.Sin(2 * math.Pi * lowFreq * float64(i) / fs) + out := pe.Process(s) + lowEnergy += out * out + } + pe.Reset() + for i := 0; i < n; i++ { + s := math.Sin(2 * math.Pi * highFreq * float64(i) / fs) + out := pe.Process(s) + highEnergy += out * out + } + + // High frequency should have more energy than low frequency + // after pre-emphasis. The input energy is the same for both. + ratio := highEnergy / lowEnergy + if ratio < 2.0 { + t.Fatalf("expected high-freq boost, energy ratio=%.2f", ratio) + } +} + +func TestPreEmphasisDeEmphasisRoundtrip(t *testing.T) { + fs := 48000.0 + pe := NewPreEmphasis(50, fs) + de := NewDeEmphasis(50, fs) + + // Run a 1kHz tone through pre+de, should approximately recover + n := 4800 + freq := 1000.0 + var maxErr float64 + + // Let filters settle for 200 samples + for i := 0; i < 200; i++ { + s := math.Sin(2 * math.Pi * freq * float64(i) / fs) + de.Process(pe.Process(s)) + } + for i := 200; i < n; i++ { + s := math.Sin(2 * math.Pi * freq * float64(i) / fs) + recovered := de.Process(pe.Process(s)) + err := math.Abs(recovered - s) + if err > maxErr { + maxErr = err + } + } + if maxErr > 0.05 { + t.Fatalf("roundtrip error too large: %.4f", maxErr) + } +} + +func TestDisabledPreEmphasis(t *testing.T) { + pe := NewPreEmphasis(0, 48000) + if pe.Process(0.5) != 0.5 { + t.Fatal("disabled filter should pass through") + } +}