From e8c5c2729b63c2a772a26d86b4d11f50ca3ba071 Mon Sep 17 00:00:00 2001 From: Jan Date: Sat, 11 Apr 2026 10:57:30 +0200 Subject: [PATCH] watermark: replace time-domain PN with STFT-domain spread-spectrum (Kirovski & Malvar 2003) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Complete redesign of the audio watermark system. The previous time-domain PN chip correlator failed over the FM air channel due to non-linear clock drift between TX and RX sample clocks, destroying bit-boundary alignment and producing ~50% BER (random). The new system embeds the watermark in STFT magnitudes (frequency domain), eliminating all sample-clock sensitivity. Based on: D. Kirovski, H.S. Malvar, "Spread-Spectrum Watermarking of Audio Signals" IEEE Trans. Signal Processing, Vol. 51, No. 4, April 2003 Architecture: Encoder (generator.go, two-pass): Audio → Stages 1-3 (LPF, drive, clip, cleanup) → Decimate ÷19 → 12 kHz mono → STFT (512-point Hann, hop=256, 21ms frames) → Magnitude × (1 ± 0.059) per PN chip (0.5 dB, multiplicative) → ISTFT → difference → ZOH ×19 → interpolation LPF@5.5kHz → add to L/R → Stages 4-6 (stereo encode, composite clip, pilot, RDS) Decoder (wmdecode, key-free): Recording → LPF@5.5kHz → decimate to 12 kHz → STFT (same parameters) → Cepstrum filter (DCT, zero first 8 coefficients → -6 dB carrier noise) → Cycle-offset search (6400 candidates, ~70s for 20-min recording) → PN correlation → 128 soft bit decisions → RS(16,8) Vandermonde decode → 64-bit fingerprint Key properties: - Multiplicative embedding: watermark scales with audio content. Loud audio masks stronger watermark, silence gets zero watermark. No gate needed (old system required silence gate to prevent audibility). - Block repetition R=5: each PN chip repeated across 5 consecutive STFT frames. Detection uses center frame only. Tolerates ±2 frames (±43ms) of timing drift without any clock recovery. - Fixed PN sequence: spreading code is public (seed "fmrtx-stft-pn-v1"), enabling blind fingerprint extraction without knowing the license key. Key identity is carried solely in the RS-encoded payload. - PCC covert channel: PN partitioned into 128 subsets (one per data bit), permuted across the watermark cycle for localized-damage resilience. Parameters: WMRate: 12000 Hz (228000/19, 192000/16 — exact both sides) FFT: 512-point, Hann window, hop=256 (50% overlap) Sub-band: bins 9-213 (200-5000 Hz), 204 frequency chips/frame Embedding: 0.5 dB (multiplicative, 0.00 dB RMS change on audio) Spreading: 204 bins × 10 groups/bit = 2040 chips/bit (33 dB gain) Block rep: R=5 (±2 frame drift tolerance) Payload: RS(16,8) → 64 bit fingerprint (SHA-256 truncated) WM cycle: 136.5 seconds Decode margin: ~19 dB over noise floor at 0.5 dB embedding Over-the-air results (PlutoSDR TX → SDRplay RX, 102.8 MHz): BER: 0/128 Erasures: 0 avg|c|: 2260 Spectrum: clean, no artifacts above 6 kHz Bug fixes included: - RS decoder: replaced Forney formula (used α^pos instead of α^(15-pos) due to polynomial convention mismatch) with Vandermonde Gaussian elimination solver. The Horner-method syndrome computation uses C(x) = c[0]x^15 + ... + c[15], mapping byte position j to polynomial power (15-j). Forney was using α^j as the error locator instead of α^(15-j), producing wrong correction magnitudes for all erasure configurations. - Generator decimation: anti-alias LPF was applied only to every 19th sample instead of all composite-rate samples. IIR filter state was updated only at 12 kHz effective rate, producing incorrect filtering. Fixed: LPF processes all 228 kHz samples, then decimate. - ZOH spectral images: zero-order hold upsample (12k→228k) created images at 12k, 24k, 36k Hz, leaking into pilot/stereo-sub/RDS bands. Fixed: separate interpolation LPF@5.5kHz on the upsample path. - Detect() single-cycle bug: iterated over groups (one frame per group) instead of all recording frames. Longer recordings did not improve SNR. Fixed: iterate all frames with modular wrapping for automatic multi-cycle averaging. New files: internal/dsp/fft.go Radix-2 Cooley-Tukey FFT/IFFT internal/watermark/stft_watermark.go STFT embedder + detector internal/watermark/stft_roundtrip_test.go Changed files: internal/watermark/watermark.go RS Vandermonde fix, accessor methods internal/offline/generator.go Two-pass architecture, STFT overlay cmd/wmdecode/main.go Complete rewrite, key-free extraction cmd/fmrtx/main.go Remove old watermark log + import Removed: Old time-domain Embedder, gate, pulse-shaping LPF, PN sequence, ChipRate/Level/RecordingRate constants, DiagnosticState --- cmd/fmrtx/main.go | 3 - cmd/wmdecode/main.go | 252 +++++---- internal/offline/generator.go | 15 +- internal/watermark/stft_roundtrip_test.go | 6 +- internal/watermark/stft_watermark.go | 36 +- internal/watermark/watermark.go | 487 +++--------------- .../watermark/watermark_roundtrip_test.go | 188 ------- 7 files changed, 222 insertions(+), 765 deletions(-) delete mode 100644 internal/watermark/watermark_roundtrip_test.go diff --git a/cmd/fmrtx/main.go b/cmd/fmrtx/main.go index 048c10b..dfdc990 100644 --- a/cmd/fmrtx/main.go +++ b/cmd/fmrtx/main.go @@ -13,7 +13,6 @@ import ( apppkg "github.com/jan/fm-rds-tx/internal/app" "github.com/jan/fm-rds-tx/internal/license" - "github.com/jan/fm-rds-tx/internal/watermark" "github.com/jan/fm-rds-tx/internal/audio" cfgpkg "github.com/jan/fm-rds-tx/internal/config" ctrlpkg "github.com/jan/fm-rds-tx/internal/control" @@ -186,8 +185,6 @@ func runTXMode(cfg cfgpkg.Config, configPath string, driver platform.SoapyDriver log.Printf("license: no valid key — evaluation jingle every %d minutes", license.JingleIntervalMinutes) } engine.SetLicenseState(licState, licenseKey) - log.Printf("watermark: STFT-domain embedding — WMRate=%d FFTSize=%d Level=%.1fdB", - watermark.WMRate, watermark.FFTSize, watermark.WMLevelDB) cfg = applyLegacyAudioFlags(cfg, audioStdin, audioRate, audioHTTP) var streamSrc *audio.StreamSource diff --git a/cmd/wmdecode/main.go b/cmd/wmdecode/main.go index 0d7cce2..ae3e487 100644 --- a/cmd/wmdecode/main.go +++ b/cmd/wmdecode/main.go @@ -1,11 +1,13 @@ // cmd/wmdecode — STFT-domain spread-spectrum watermark decoder. // -// Decodes watermark from FM broadcast recordings following -// Kirovski & Malvar (IEEE TSP 2003) architecture. +// Extracts the embedded fingerprint from an FM broadcast recording +// without needing to know the license key. Optionally verifies +// against supplied keys. // // Usage: // -// wmdecode [key ...] +// wmdecode Extract fingerprint +// wmdecode [key ...] Extract + verify keys package main import ( @@ -55,9 +57,8 @@ func main() { for i := 0; i < nDown; i++ { down[i] = filtered[i*decimFactor] } - fmt.Printf("Downsampled: %d samples, %.1fs\n", nDown, float64(nDown)/wmRate) - // Step 2: Compute ALL STFT frames with cepstrum filtering + // Step 2: STFT + cepstrum filtering fftSize := watermark.FFTSize hop := watermark.FFTHop nFrames := (nDown - fftSize) / hop @@ -68,7 +69,8 @@ func main() { var window [watermark.FFTSize]float64 dsp.HannWindow(window[:]) - fmt.Printf("STFT: %d frames (%d-point, hop=%d)\n", nFrames, fftSize, hop) + fmt.Printf("STFT: %d frames (%.1fs required, %.1fs available)\n", + nFrames, float64(watermark.SamplesPerWM)/wmRate, float64(nDown)/wmRate) type stftMag [watermark.FFTSize / 2]float64 frameMags := make([]stftMag, nFrames) @@ -89,151 +91,137 @@ func main() { cepstrumFilter(frameMags[f][:], 8) } - // Step 3: For each key, search cycle offset + rep offset - keys := os.Args[2:] - if len(keys) == 0 { - fmt.Println("No keys supplied.") - os.Exit(1) - } - - for _, key := range keys { - fmt.Printf("\nKey: %q\n", key) - det := watermark.NewSTFTDetector(key) - - totalGroups := watermark.TotalGroups - timeRep := watermark.TimeRep - framesPerWM := watermark.FramesPerWM - numBins := watermark.NumBins - binLow := watermark.BinLow - centerRep := timeRep / 2 - - bestMetric := -1.0 - var bestCorrs [watermark.PayloadBits]float64 - bestCycleOff := 0 - bestRepOff := 0 - - nCandidates := 0 - for cycleOff := 0; cycleOff < framesPerWM; cycleOff += timeRep { - for repOff := 0; repOff < timeRep; repOff++ { - var testCorrs [watermark.PayloadBits]float64 - - for f := 0; f < nFrames; f++ { - wmFrame := ((f - cycleOff - repOff) % framesPerWM + framesPerWM) % framesPerWM - if wmFrame%timeRep != centerRep { - continue - } - g := wmFrame / timeRep - if g >= totalGroups { - continue - } - - var corr float64 - for b := 0; b < numBins; b++ { - corr += frameMags[f][binLow+b] * float64(det.PNChipAt(g, b)) - } - testCorrs[det.GroupBit(g)] += corr + // Step 3: Cycle-offset search (key-independent — fixed PN) + det := watermark.NewSTFTDetector() + + totalGroups := watermark.TotalGroups + timeRep := watermark.TimeRep + framesPerWM := watermark.FramesPerWM + numBins := watermark.NumBins + binLow := watermark.BinLow + centerRep := timeRep / 2 + + bestMetric := -1.0 + var bestCorrs [watermark.PayloadBits]float64 + bestCycleOff := 0 + bestRepOff := 0 + + for cycleOff := 0; cycleOff < framesPerWM; cycleOff += timeRep { + for repOff := 0; repOff < timeRep; repOff++ { + var testCorrs [watermark.PayloadBits]float64 + for f := 0; f < nFrames; f++ { + wmFrame := ((f - cycleOff - repOff) % framesPerWM + framesPerWM) % framesPerWM + if wmFrame%timeRep != centerRep { + continue } - - var metric float64 - for _, c := range testCorrs { - metric += c * c + g := wmFrame / timeRep + if g >= totalGroups { + continue } - - if metric > bestMetric { - bestMetric = metric - bestCorrs = testCorrs - bestCycleOff = cycleOff - bestRepOff = repOff + var corr float64 + for b := 0; b < numBins; b++ { + corr += frameMags[f][binLow+b] * float64(det.PNChipAt(g, b)) } - nCandidates++ + testCorrs[det.GroupBit(g)] += corr } - } - - fmt.Printf("Searched %d candidates in %v\n", nCandidates, time.Since(t0).Round(time.Millisecond)) - fmt.Printf("Best: cycleOff=%d, repOff=%d, metric=%.0f\n", bestCycleOff, bestRepOff, bestMetric) - - var sumAbs float64 - for _, c := range bestCorrs { - sumAbs += math.Abs(c) - } - fmt.Printf("Corrs: avg|c|=%.1f\n", sumAbs/128) - - // BER diagnostic against known key - knownPayload := watermark.KeyToPayload(key) - knownCW := watermark.RSEncode(knownPayload) - var knownBits [watermark.PayloadBits]int - for i := 0; i < watermark.PayloadBits; i++ { - knownBits[i] = int((knownCW[i/8] >> uint(7-(i%8))) & 1) - } - nerr := 0 - for i := 0; i < watermark.PayloadBits; i++ { - hard := 0 - if bestCorrs[i] < 0 { - hard = 1 + var metric float64 + for _, c := range testCorrs { + metric += c * c } - if hard != knownBits[i] { - nerr++ + if metric > bestMetric { + bestMetric = metric + bestCorrs = testCorrs + bestCycleOff = cycleOff + bestRepOff = repOff } } - fmt.Printf("BER: %d/128 (%.1f%%)\n", nerr, 100*float64(nerr)/128) - - // Show recv vs expected - var recv [watermark.RsTotalBytes]byte - confs := make([]float64, watermark.PayloadBits) - for i := 0; i < watermark.PayloadBits; i++ { - confs[i] = math.Abs(bestCorrs[i]) - if bestCorrs[i] < 0 { - recv[i/8] |= 1 << uint(7-(i%8)) - } + } + + var sumAbs float64 + for _, c := range bestCorrs { + sumAbs += math.Abs(c) + } + fmt.Printf("Sync: cycleOff=%d repOff=%d (%.1fs into recording)\n", + bestCycleOff, bestRepOff, float64(bestCycleOff*hop)/wmRate) + fmt.Printf("Corrs: avg|c|=%.1f\n", sumAbs/128) + + // Step 4: RS decode — extract fingerprint + var recv [watermark.RsTotalBytes]byte + confs := make([]float64, watermark.PayloadBits) + for i := 0; i < watermark.PayloadBits; i++ { + confs[i] = math.Abs(bestCorrs[i]) + if bestCorrs[i] < 0 { + recv[i/8] |= 1 << uint(7-(i%8)) } - fmt.Printf("recv: %x\nwant: %x\n", recv, knownCW) - - // Confidence-based erasure (MIN bit confidence per byte) - type bc struct{ idx int; conf float64 } - byteConfs := make([]bc, watermark.RsTotalBytes) - for b := 0; b < watermark.RsTotalBytes; b++ { - minC := confs[b*8] - for bit := 1; bit < 8; bit++ { - if confs[b*8+bit] < minC { - minC = confs[b*8+bit] - } + } + + type bc struct{ idx int; conf float64 } + byteConfs := make([]bc, watermark.RsTotalBytes) + for b := 0; b < watermark.RsTotalBytes; b++ { + minC := confs[b*8] + for bit := 1; bit < 8; bit++ { + if confs[b*8+bit] < minC { + minC = confs[b*8+bit] } - byteConfs[b] = bc{b, minC} } - sort.Slice(byteConfs, func(a, b int) bool { return byteConfs[a].conf < byteConfs[b].conf }) - - decoded := false - for nErase := 0; nErase <= watermark.RsCheckBytes; nErase++ { - if nErase == 0 { - p, ok := watermark.RSDecode(recv, nil) - if ok && watermark.KeyMatchesPayload(key, p) { - fmt.Printf(" ✓ MATCH (0 erasures), payload=%x\n", p) - decoded = true - break - } - continue - } - erasePos := make([]int, nErase) - for i := 0; i < nErase; i++ { - erasePos[i] = byteConfs[i].idx - } - sort.Ints(erasePos) - p, ok := watermark.RSDecode(recv, erasePos) - if ok && watermark.KeyMatchesPayload(key, p) { - fmt.Printf(" ✓ MATCH (%d erasures), payload=%x\n", nErase, p) + byteConfs[b] = bc{b, minC} + } + sort.Slice(byteConfs, func(a, b int) bool { return byteConfs[a].conf < byteConfs[b].conf }) + + var payload [watermark.RsDataBytes]byte + decoded := false + nErasures := 0 + for nErase := 0; nErase <= watermark.RsCheckBytes; nErase++ { + if nErase == 0 { + p, ok := watermark.RSDecode(recv, nil) + if ok { + payload = p decoded = true break } + continue + } + erasePos := make([]int, nErase) + for i := 0; i < nErase; i++ { + erasePos[i] = byteConfs[i].idx } + sort.Ints(erasePos) + p, ok := watermark.RSDecode(recv, erasePos) + if ok { + payload = p + nErasures = nErase + decoded = true + break + } + } - if !decoded { - fmt.Println(" ✗ NOT FOUND") + if !decoded { + fmt.Printf("\nWatermark: NOT DETECTED\n") + fmt.Printf("Done in %v\n", time.Since(t0).Round(time.Millisecond)) + os.Exit(1) + } + + fmt.Printf("\nWatermark: DETECTED (%d erasures)\n", nErasures) + fmt.Printf("Fingerprint: %x\n", payload) + + // Optional: verify against supplied keys + keys := os.Args[2:] + if len(keys) > 0 { + fmt.Println() + for _, key := range keys { + if watermark.KeyMatchesPayload(key, payload) { + fmt.Printf(" ✓ MATCH: %q\n", key) + } else { + fmt.Printf(" ✗ %q\n", key) + } } } fmt.Printf("\nDone in %v\n", time.Since(t0).Round(time.Millisecond)) } +// --- Cepstrum filter --- + func cepstrumFilter(magDB []float64, nCeps int) { n := len(magDB) if n < nCeps*2 { @@ -263,6 +251,8 @@ func cepstrumFilter(magDB []float64, nCeps int) { } } +// --- LPF --- + type biquad struct{ b0, b1, b2, a1, a2 float64 } type iirCoeffs []biquad @@ -299,6 +289,8 @@ func applyIIR(samples []float64, coeffs iirCoeffs) []float64 { return out } +// --- WAV reader --- + func rmsLevel(s []float64) float64 { var acc float64 for _, v := range s { diff --git a/internal/offline/generator.go b/internal/offline/generator.go index ccce1d0..4a0bb20 100644 --- a/internal/offline/generator.go +++ b/internal/offline/generator.go @@ -135,6 +135,7 @@ type Generator struct { // Watermark: STFT-domain spread-spectrum (Kirovski & Malvar 2003). stftEmbedder *watermark.STFTEmbedder wmDecimLPF *dsp.FilterChain // anti-alias LPF for 228k→12k decimation + wmInterpLPF *dsp.FilterChain // image-rejection LPF for 12k→228k upsample } func NewGenerator(cfg cfgpkg.Config) *Generator { @@ -284,6 +285,7 @@ func (g *Generator) init() { // Nyquist at 12 kHz = 6 kHz. Cut at 5.5 kHz with margin. if g.stftEmbedder != nil { g.wmDecimLPF = dsp.NewLPF4(5500, g.sampleRate) + g.wmInterpLPF = dsp.NewLPF4(5500, g.sampleRate) // separate instance for upsample } g.initialized = true @@ -447,13 +449,19 @@ func (g *Generator) GenerateFrame(duration time.Duration) *output.CompositeFrame // STFT embed at 12 kHz embedded := g.stftEmbedder.ProcessBlock(mono12k) - // Extract watermark signal (difference) and upsample via ZOH + // Extract watermark signal (difference) and upsample via ZOH + LPF. + // ZOH creates spectral images at 12k, 24k, 36k... Hz. + // The interpolation LPF removes these, keeping only 0-5.5 kHz. + // Without this, the images leak into pilot (19k) and stereo sub (38k). for i := 0; i < samples; i++ { wmIdx := i / decimFactor if wmIdx >= nDown { wmIdx = nDown - 1 } wmSig := embedded[wmIdx] - mono12k[wmIdx] + if g.wmInterpLPF != nil { + wmSig = g.wmInterpLPF.Process(wmSig) + } lBuf[i] += wmSig rBuf[i] += wmSig } @@ -523,11 +531,6 @@ func (g *Generator) GenerateFrame(duration time.Duration) *output.CompositeFrame g.bs412.ProcessChunk(bs412PowerAccum / float64(samples)) } - // STFT watermark diagnostic - if g.stftEmbedder != nil && g.frameSeq%100 == 1 { - log.Printf("watermark stft: frame=%d, active", g.frameSeq) - } - return frame } diff --git a/internal/watermark/stft_roundtrip_test.go b/internal/watermark/stft_roundtrip_test.go index c91b081..f23e60d 100644 --- a/internal/watermark/stft_roundtrip_test.go +++ b/internal/watermark/stft_roundtrip_test.go @@ -36,8 +36,9 @@ func TestSTFTRoundTrip(t *testing.T) { t.Logf("RMS change: %.2f dB", 20*math.Log10(rmsOut/rmsIn)) // Detect watermark - detector := NewSTFTDetector(key) - corrs, offset := detector.Detect(watermarked) + detector := NewSTFTDetector() + // Skip embedder OLA startup (first FFTSize samples are pass-through, not watermarked) + corrs, offset := detector.Detect(watermarked[FFTSize:]) t.Logf("Detection offset: %d", offset) @@ -75,6 +76,7 @@ func TestSTFTRoundTrip(t *testing.T) { } if hard != expectedBits[i] { nerr++ + t.Logf(" ERR bit %d: corr=%+.2f, expected=%d got=%d", i, corrs[i], expectedBits[i], hard) } } t.Logf("BER: %d/%d (%.1f%%)", nerr, payloadBits, 100*float64(nerr)/float64(payloadBits)) diff --git a/internal/watermark/stft_watermark.go b/internal/watermark/stft_watermark.go index 9548d15..6595d5d 100644 --- a/internal/watermark/stft_watermark.go +++ b/internal/watermark/stft_watermark.go @@ -32,7 +32,7 @@ const ( NumBins = BinHigh - BinLow // 204 frequency chips per STFT frame TimeRep = 5 // block repetition factor (±2 frame drift tolerance) GroupsPerBit = 10 // time groups per data bit - WMLevelDB = 1.5 // embedding level (dB) + WMLevelDB = 0.5 // embedding level (dB) — inaudible, 20 dB margin for decode TotalGroups = GroupsPerBit * payloadBits // 10 × 128 = 1280 FramesPerWM = TotalGroups * TimeRep // 1280 × 5 = 6400 @@ -63,7 +63,7 @@ type STFTEmbedder struct { frameIdx int // STFT frame counter (wraps at FramesPerWM) primed bool // true after first full frame - // Level in linear scale: 10^(WMLevelDB/20) - 1 ≈ 0.189 for 1.5 dB + // Level in linear scale: 10^(WMLevelDB/20) - 1 ≈ 0.059 for 0.5 dB levelLinear float64 } @@ -89,7 +89,7 @@ func NewSTFTEmbedder(key string) *STFTEmbedder { } // Generate PN chips from key-seeded PRNG - seed := sha256.Sum256(append([]byte("stft-pn-"), key...)) + seed := sha256.Sum256([]byte("fmrtx-stft-pn-v1")) prng := newPRNG(seed[:]) for g := 0; g < TotalGroups; g++ { for b := 0; b < NumBins; b++ { @@ -107,7 +107,7 @@ func NewSTFTEmbedder(key string) *STFTEmbedder { e.groupToBit[g] = g % payloadBits } // Permute within each bit's groups using key-seeded PRNG - permSeed := sha256.Sum256(append([]byte("stft-perm-"), key...)) + permSeed := sha256.Sum256([]byte("fmrtx-stft-perm-v1")) permRNG := newPRNG(permSeed[:]) for i := TotalGroups - 1; i > 0; i-- { j := permRNG.next() % uint32(i+1) @@ -232,12 +232,13 @@ type STFTDetector struct { groupToBit [TotalGroups]int } -// NewSTFTDetector creates a detector matching the given key's PN sequence. -func NewSTFTDetector(key string) *STFTDetector { +// NewSTFTDetector creates a detector. No key needed — the PN sequence is +// public (fixed). The detector extracts the payload blindly. +func NewSTFTDetector() *STFTDetector { d := &STFTDetector{} // Same PN generation as embedder - seed := sha256.Sum256(append([]byte("stft-pn-"), key...)) + seed := sha256.Sum256([]byte("fmrtx-stft-pn-v1")) prng := newPRNG(seed[:]) for g := 0; g < TotalGroups; g++ { for b := 0; b < NumBins; b++ { @@ -253,7 +254,7 @@ func NewSTFTDetector(key string) *STFTDetector { for g := 0; g < TotalGroups; g++ { d.groupToBit[g] = g % payloadBits } - permSeed := sha256.Sum256(append([]byte("stft-perm-"), key...)) + permSeed := sha256.Sum256([]byte("fmrtx-stft-perm-v1")) permRNG := newPRNG(permSeed[:]) for i := TotalGroups - 1; i > 0; i-- { j := permRNG.next() % uint32(i+1) @@ -317,21 +318,24 @@ func (d *STFTDetector) Detect(audio []float64) (corrs [payloadBits]float64, best for startOffset := 0; startOffset < TimeRep; startOffset++ { var testCorrs [payloadBits]float64 - // For each group, use the CENTER frame of the repetition block - for g := 0; g < TotalGroups; g++ { - bitIdx := d.groupToBit[g] - frameInWM := g*TimeRep + startOffset + TimeRep/2 - if frameInWM >= nFrames { + // Iterate over ALL recording frames — multiple WM cycles accumulate + // automatically via modular wrapping. This gives √N_cycles SNR gain. + for f := 0; f < nFrames; f++ { + wmFrame := ((f - startOffset) % FramesPerWM + FramesPerWM) % FramesPerWM + if wmFrame%TimeRep != TimeRep/2 { + continue // not center of repetition block + } + g := wmFrame / TimeRep + if g >= TotalGroups { continue } - // Correlate this frame's magnitudes with the PN chips var corr float64 for b := 0; b < NumBins; b++ { bin := BinLow + b - corr += frames[frameInWM].magDB[bin] * float64(d.pnChips[g][b]) + corr += frames[f].magDB[bin] * float64(d.pnChips[g][b]) } - testCorrs[bitIdx] += corr + testCorrs[d.groupToBit[g]] += corr } // Detection metric: sum of squared partial correlations (chi-squared) diff --git a/internal/watermark/watermark.go b/internal/watermark/watermark.go index 571997c..1e2feb7 100644 --- a/internal/watermark/watermark.go +++ b/internal/watermark/watermark.go @@ -1,223 +1,37 @@ -// Package watermark implements spread-spectrum audio watermarking for fm-rds-tx. +// Package watermark implements STFT-domain spread-spectrum audio watermarking +// for fm-rds-tx, based on Kirovski & Malvar (IEEE TSP 2003). // -// # Design +// The watermark is embedded in STFT magnitude (dB scale) — a multiplicative +// modification that naturally scales with audio content (psychoacoustic masking). +// Block repetition coding provides drift tolerance. Reed-Solomon RS(16,8) +// provides error correction for the 128-bit payload. // -// The watermark is injected into the audio L/R signal after all audio -// processing (LPF, clip, limiter) but before stereo encoding, so it -// survives FM broadcast, receiver demodulation, de-emphasis, and moderate -// EQ intact. The PN chip rate is limited to 12 kHz so all watermark energy -// sits within 0–6 kHz — the band that survives every stage of the FM chain: -// -// Embedder → Stereo Encode → Composite Clip → Pilot/RDS → FM Mod -// → FM Demod → De-emphasis → Stereo Decode → Audio Out → Recording -// -// The payload is Reed-Solomon encoded for robust recovery even when -// individual bits have high error rates due to noise and audio masking. -// -// # Parameters -// -// - PN sequence: 2048-chip LFSR-13 (seed 0x1ACE) -// - Payload: 8 bytes (SHA-256[:8] of key) → RS(16,8) → 16 bytes → 128 bits -// - Chip clock: 12 kHz → PN bandwidth 0–6 kHz (survives de-emphasis) -// - Frame period: ~21.8 s at 228 kHz composite (repeats ~2.7×/min) -// - Injection: -48 dBFS on audio L+R after all processing, before stereo encode -// - Spreading gain: 33 dB. RS erasure corrects up to 8 of 16 byte symbols. -// -// # Recovery (cmd/wmdecode) -// -// 1. Record FM receiver audio output as mono WAV (any sample rate ≥ 12 kHz). -// 2. Phase search: energy-based coarse/fine search for chip alignment. -// 3. Extract 128 bit correlations at found phase, averaged over all frames. -// 4. Frame sync: try all 128 cyclic rotations of the bit sequence, -// RS-decode each; the rotation that succeeds gives the frame alignment. -// 5. Byte-level erasure + soft-decision bit-flipping for error correction. -// 6. RS erasure-decode → 8 payload bytes → compare against known keys. +// Encoder: internal/watermark/stft_watermark.go (STFTEmbedder) +// Decoder: cmd/wmdecode (STFTDetector) package watermark import ( "crypto/sha256" + "fmt" ) +// RS code parameters. const ( - // pnChips is the spreading factor — PN chips per data bit. - // Spreading gain = 10·log10(2048) = 33.1 dB. - pnChips = 2048 - - // rsDataBytes is the number of payload bytes before RS encoding. - rsDataBytes = 8 - - // rsCheckBytes is the number of RS parity bytes. With 8 check bytes the - // code corrects up to 4 errors or up to 8 erasures per 16-byte codeword. - rsCheckBytes = 8 - - // rsTotalBytes is the full RS codeword length. + rsDataBytes = 8 // payload bytes + rsCheckBytes = 8 // parity bytes rsTotalBytes = rsDataBytes + rsCheckBytes // 16 - - // payloadBits is the total number of BPSK bits per watermark frame. - payloadBits = rsTotalBytes * 8 // 128 - - // Level is the audio injection amplitude per channel (-48 dBFS). - // At typical audio levels this is completely inaudible. - Level = 0.040 - - // CompositeRate is the sample rate at which the watermark is embedded. - CompositeRate = 228000 + payloadBits = rsTotalBytes * 8 // 128 ) -// ChipRate is the effective PN chip clock rate (Hz). Determines the spectral -// bandwidth of the watermark: Nyquist = ChipRate/2 = 6 kHz. This ensures -// all PN energy is within the audio band that survives de-emphasis (50/75 µs), -// receiver LPFs, audio codecs, speaker EQ, and even acoustic recording. -// -// At CompositeRate (228 kHz), each chip spans 228000/12000 = 19 samples. -// At any recording rate R, each chip spans R/12000 samples. -const ChipRate = 12000 - -// RecordingRate is the canonical recording rate for test WAV output (wmtest). -// Not used for chip stepping — ChipRate controls that. -const RecordingRate = 48000 - -// Embedder continuously embeds a watermark into audio L/R samples. -// Not thread-safe: call NextSample from the single DSP goroutine only. -type Embedder struct { - codeword [rsTotalBytes]byte // RS-encoded payload, 16 bytes - chipIdx int // chip position within current bit (0..pnChips-1) - bitIdx int // current bit in codeword (0..127) - symbol int8 // BPSK symbol for current bit: +1 or -1 - accum int // Bresenham accumulator for chip-rate stepping - - // Audio-level gate: mutes watermark during silence to prevent audibility. - gateGain float64 // smooth ramp 0.0 (muted) → 1.0 (open) - gateThreshold float64 // audio level below which gate closes - gateRampUp float64 // per-sample increment when opening (~5ms) - gateRampDown float64 // per-sample decrement when closing (~5ms) - gateEnabled bool -} - -// NewEmbedder creates an Embedder for the given license key. -// The key's SHA-256 hash (first 8 bytes) is RS-encoded and embedded. -// An empty key embeds a null payload (still watermarks, just anonymous). -func NewEmbedder(key string) *Embedder { - var data [rsDataBytes]byte - if key != "" { - h := sha256.Sum256([]byte(key)) - copy(data[:], h[:rsDataBytes]) - } - e := &Embedder{gateGain: 1.0} - e.codeword = rsEncode(data) - e.loadSymbol() - return e -} - -// NextSample returns the watermark amplitude for one composite sample. -// Add this value to both audio.Frame.L and audio.Frame.R before stereo encoding. -// -// The chip index advances using Bresenham stepping at ChipRate/CompositeRate, -// so each chip occupies exactly CompositeRate/ChipRate composite samples on -// average (~19 samples at 228 kHz). The PN signal bandwidth is 0–6 kHz. -func (e *Embedder) NextSample() float64 { - chip := float64(pnSequence[e.chipIdx]) - sample := Level * float64(e.symbol) * chip * e.gateGain - - // Bresenham: advance chip once per ChipRate/CompositeRate composite samples. - e.accum += ChipRate - if e.accum >= CompositeRate { - e.accum -= CompositeRate - e.chipIdx++ - if e.chipIdx >= pnChips { - e.chipIdx = 0 - e.bitIdx = (e.bitIdx + 1) % payloadBits - e.loadSymbol() - } - } - return sample -} - -// loadSymbol sets e.symbol from the current bit in the codeword (MSB first). -func (e *Embedder) loadSymbol() { - byteIdx := e.bitIdx / 8 - bitPos := uint(7 - (e.bitIdx % 8)) - if (e.codeword[byteIdx]>>bitPos)&1 == 0 { - e.symbol = 1 - } else { - e.symbol = -1 - } -} - -// PayloadHex returns the RS-encoded codeword as hex for logging. -func (e *Embedder) PayloadHex() string { - const hx = "0123456789abcdef" - out := make([]byte, rsTotalBytes*2) - for i, b := range e.codeword { - out[i*2] = hx[b>>4] - out[i*2+1] = hx[b&0xf] - } - return string(out) -} - -// EnableGate activates audio-level gating with asymmetric ramp times. -// threshold is the linear audio amplitude below which the watermark is muted -// (e.g. 0.01 ≈ -40 dBFS). compositeRate is needed to compute ramp speed. -// Attack (open) is fast (5ms) so the watermark starts immediately with audio. -// Release (close) is slow (200ms) to keep the watermark running through normal -// inter-word and inter-phrase gaps — only extended silence mutes. -func (e *Embedder) EnableGate(threshold, compositeRate float64) { - attackSamples := compositeRate * 0.005 // 5ms open - releaseSamples := compositeRate * 0.200 // 200ms close - if attackSamples < 1 { - attackSamples = 1 - } - if releaseSamples < 1 { - releaseSamples = 1 - } - e.gateThreshold = threshold - e.gateRampUp = 1.0 / attackSamples - e.gateRampDown = 1.0 / releaseSamples - e.gateEnabled = true -} - -// SetAudioLevel updates the gate state based on the current audio amplitude. -// Call once per sample before NextSample. absLevel should be the absolute -// mono audio level (pre- or post-pre-emphasis, either works). -func (e *Embedder) SetAudioLevel(absLevel float64) { - if !e.gateEnabled { - return - } - if absLevel > e.gateThreshold { - e.gateGain += e.gateRampUp - if e.gateGain > 1.0 { - e.gateGain = 1.0 - } - } else { - e.gateGain -= e.gateRampDown - if e.gateGain < 0.0 { - e.gateGain = 0.0 - } - } -} - -// DiagnosticState returns internal state for debugging. -type DiagnosticInfo struct { - GateGain float64 - GateEnabled bool - ChipIdx int - BitIdx int - Symbol int8 -} - -// DiagnosticState returns a snapshot of the embedder's internal state. -func (e *Embedder) DiagnosticState() DiagnosticInfo { - return DiagnosticInfo{ - GateGain: e.gateGain, - GateEnabled: e.gateEnabled, - ChipIdx: e.chipIdx, - BitIdx: e.bitIdx, - Symbol: e.symbol, - } -} +// Exported constants for tools. +const ( + PayloadBits = payloadBits + RsDataBytes = rsDataBytes + RsTotalBytes = rsTotalBytes + RsCheckBytes = rsCheckBytes +) -// --- RS(16,8) over GF(2^8) — GF poly 0x11d, fcr=0, generator=2 --- -// These routines are used by the embedder (encode) and the recovery tool (decode). +// --- GF(2^8) arithmetic — primitive polynomial 0x11d --- func gfMul(a, b byte) byte { if a == 0 || b == 0 { @@ -234,17 +48,21 @@ func gfInv(a byte) byte { } func gfPow(a byte, n int) byte { - if a == 0 { - return 0 + if n == 0 { + return 1 } - return gfExp[(int(gfLog[a])*n)%255] + r := a + for i := 1; i < n; i++ { + r = gfMul(r, a) + } + return r } -// rsEncode encodes 8 data bytes into a 16-byte RS codeword. +// --- RS(16,8) encoding --- + func rsEncode(data [rsDataBytes]byte) [rsTotalBytes]byte { - var work [rsTotalBytes]byte - copy(work[:rsDataBytes], data[:]) - // Polynomial long division by the generator polynomial. + work := make([]byte, rsTotalBytes) + copy(work, data[:]) for i := 0; i < rsDataBytes; i++ { fb := work[i] if fb != 0 { @@ -253,20 +71,29 @@ func rsEncode(data [rsDataBytes]byte) [rsTotalBytes]byte { } } } - var cw [rsTotalBytes]byte - copy(cw[:rsDataBytes], data[:]) - copy(cw[rsDataBytes:], work[rsDataBytes:]) - return cw + var out [rsTotalBytes]byte + copy(out[:rsDataBytes], data[:]) + copy(out[rsDataBytes:], work[rsDataBytes:]) + return out +} + +// RSEncode encodes 8 data bytes into a 16-byte RS codeword (exported). +func RSEncode(data [rsDataBytes]byte) [rsTotalBytes]byte { + return rsEncode(data) } +// --- RS(16,8) decoding (Vandermonde solver) --- + // RSDecode recovers 8 data bytes from a (possibly corrupted) 16-byte codeword. -// erasurePositions lists the byte indices (0..15) of symbols with low confidence -// that should be treated as erasures. Up to 8 erasures can be corrected. -// Returns (data, true) on success, (zero, false) on decoding failure. +// erasurePositions lists the byte indices (0-15) of known-erased bytes. +// Returns the 8 payload bytes and true if decoding succeeded. +// +// The polynomial convention is C(x) = c[0]x^15 + c[1]x^14 + … + c[15], +// so byte position j maps to polynomial power (15-j). The locator for +// position j is α^(15-j). Decoding uses Vandermonde/Gaussian elimination +// instead of Forney's formula to avoid position-mapping bugs. func RSDecode(recv [rsTotalBytes]byte, erasurePositions []int) ([rsDataBytes]byte, bool) { - // Step 1: compute syndromes S[i] = C(α^i) via Horner's method. - // The polynomial convention is C(x) = c[0]x^15 + c[1]x^14 + … + c[15], - // so byte position j contributes c[j]·(α^i)^(15-j) to S[i]. + // Syndromes S[i] = C(α^i) via Horner's method. var S [rsCheckBytes]byte for i := 0; i < rsCheckBytes; i++ { var acc byte @@ -276,7 +103,7 @@ func RSDecode(recv [rsTotalBytes]byte, erasurePositions []int) ([rsDataBytes]byt S[i] = acc } - // Step 2: all syndromes zero → valid codeword. + // All syndromes zero → valid codeword. hasErr := false for _, s := range S { if s != 0 { @@ -294,23 +121,13 @@ func RSDecode(recv [rsTotalBytes]byte, erasurePositions []int) ([rsDataBytes]byt return [rsDataBytes]byte{}, false } - // Step 3: Solve for error magnitudes via Vandermonde system. - // - // Because C(x) = c[0]x^15 + … + c[15], byte position j maps to - // polynomial power (15-j). The "locator" for position j is - // X_j = α^(15-j). The syndrome equation becomes: - // - // S[i] = Σ_k e_k · X_k^i - // - // This is a linear system (Vandermonde) in the unknowns e_k. - // Solve by Gaussian elimination in GF(2^8). - - // Build locators + // Locators: X[k] = α^(15-pos) for byte position pos. X := make([]byte, ne) for k, pos := range erasurePositions { X[k] = gfPow(2, rsTotalBytes-1-pos) } + // Vandermonde system: S[i] = Σ_k e_k · X[k]^i // Augmented matrix [V | S], ne × (ne+1) mat := make([][]byte, ne) for i := range mat { @@ -321,7 +138,7 @@ func RSDecode(recv [rsTotalBytes]byte, erasurePositions []int) ([rsDataBytes]byt mat[i][ne] = S[i] } - // Gaussian elimination with partial pivoting + // Gaussian elimination for col := 0; col < ne; col++ { pivot := -1 for row := col; row < ne; row++ { @@ -331,7 +148,7 @@ func RSDecode(recv [rsTotalBytes]byte, erasurePositions []int) ([rsDataBytes]byt } } if pivot < 0 { - return [rsDataBytes]byte{}, false // singular + return [rsDataBytes]byte{}, false } mat[col], mat[pivot] = mat[pivot], mat[col] inv := gfInv(mat[col][col]) @@ -358,7 +175,7 @@ func RSDecode(recv [rsTotalBytes]byte, erasurePositions []int) ([rsDataBytes]byt result[erasurePositions[k]] ^= mat[k][ne] } - // Verify syndromes after correction + // Verify for i := 0; i < rsCheckBytes; i++ { var acc byte for _, c := range result { @@ -374,6 +191,8 @@ func RSDecode(recv [rsTotalBytes]byte, erasurePositions []int) ([rsDataBytes]byt return out, true } +// --- Key utilities --- + // KeyMatchesPayload returns true if SHA-256(key)[:8] matches payload. func KeyMatchesPayload(key string, payload [rsDataBytes]byte) bool { h := sha256.Sum256([]byte(key)) @@ -390,11 +209,15 @@ func KeyToPayload(key string) [rsDataBytes]byte { return data } -// RSEncode encodes 8 data bytes into a 16-byte RS codeword. -func RSEncode(data [rsDataBytes]byte) [rsTotalBytes]byte { - return rsEncode(data) +// PayloadHex returns the hex string of the RS-encoded payload for a key. +func PayloadHex(key string) string { + data := KeyToPayload(key) + cw := rsEncode(data) + return fmt.Sprintf("%x", cw) } +// --- STFT detector accessors --- + // PNChipAt returns the PN chip value at group g, bin b. func (d *STFTDetector) PNChipAt(g, b int) int8 { return d.pnChips[g][b] @@ -405,184 +228,8 @@ func (d *STFTDetector) GroupBit(g int) int { return d.groupToBit[g] } -// Constants exported for the recovery tool and legacy tools. -const ( - PnChips = pnChips - PayloadBits = payloadBits - RsDataBytes = rsDataBytes - RsTotalBytes = rsTotalBytes - RsCheckBytes = rsCheckBytes -) - -// PnSequenceAt returns the PN chip value (+1.0 or -1.0) at the given index. -// Used by the decoder for rate-compensated correlation. -func PnSequenceAt(chipIdx int) float64 { - return float64(pnSequence[chipIdx%pnChips]) -} - -// PNSequence exposes the raw PN chip values for the chip-rate decoder. -var PNSequence = &pnSequence - -// CorrelateAt returns the correlation of received samples at the given bit -// position. recRate is the WAV sample rate. -// -// At any recording rate, chips are mapped via ChipRate: each chip spans -// recRate/ChipRate samples. At 48 kHz: 4 samples/chip, 8192 samples/bit. -// At 192 kHz: 16 samples/chip, 32768 samples/bit. -func CorrelateAt(samples []float64, bitStart int, recRate float64) float64 { - samplesPerBit := int(float64(pnChips) * recRate / float64(ChipRate)) - if samplesPerBit < 1 { - samplesPerBit = 1 - } - n := samplesPerBit - if bitStart+n > len(samples) { - n = len(samples) - bitStart - } - var acc float64 - for i := 0; i < n; i++ { - chipIdx := int(float64(i)*float64(ChipRate)/recRate) % pnChips - acc += samples[bitStart+i] * float64(pnSequence[chipIdx]) - } - return acc -} - - -// pnSequence is the 2048-chip LFSR-13 spreading code (seed 0x1ACE). -var pnSequence = [pnChips]int8{ - 1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, 1, -1, -1, - -1, -1, -1, 1, -1, 1, 1, -1, 1, 1, -1, 1, 1, 1, -1, 1, - -1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, - 1, -1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, - -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1, 1, - -1, 1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, - 1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, - -1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, - -1, 1, -1, 1, -1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, - -1, -1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, 1, 1, -1, - 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, 1, - -1, -1, 1, -1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, 1, -1, - -1, -1, -1, 1, 1, -1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, - 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 1, -1, -1, - -1, 1, -1, 1, -1, 1, 1, -1, 1, 1, -1, 1, 1, 1, 1, -1, - -1, -1, -1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, - 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, - -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, 1, 1, 1, 1, - -1, 1, -1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, - 1, 1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, - -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, 1, -1, 1, - -1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, - -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, 1, 1, 1, 1, - -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, - 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1, -1, -1, - -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, - 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, -1, -1, 1, - 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, - 1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, - 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, - -1, 1, -1, 1, -1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, - 1, 1, -1, 1, 1, 1, -1, 1, -1, -1, -1, -1, 1, 1, -1, 1, - -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, - -1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, - -1, 1, 1, 1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, - 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, - -1, 1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, - 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, -1, - 1, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, 1, - -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, - -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, - 1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, - 1, -1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, - 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, - 1, -1, -1, -1, 1, 1, 1, -1, 1, -1, 1, 1, -1, -1, -1, 1, - -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, -1, -1, -1, -1, 1, - 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, -1, 1, -1, 1, - 1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, - -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, -1, - 1, -1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, - -1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, - 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, - 1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, - 1, 1, -1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, -1, -1, 1, - 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, 1, 1, 1, - 1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, - 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, - 1, -1, -1, 1, -1, -1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, - -1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, 1, 1, - -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, -1, 1, - 1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, -1, - -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, - 1, 1, 1, 1, -1, -1, -1, 1, -1, 1, 1, 1, 1, 1, 1, 1, - 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, - -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, - -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, 1, -1, - 1, -1, -1, -1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, - 1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, - -1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, 1, -1, - -1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, 1, -1, 1, -1, - 1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, - -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, - -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, -1, -1, - 1, -1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, - -1, 1, 1, 1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, 1, - 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, - -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, -1, 1, - -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, - -1, 1, 1, -1, 1, -1, 1, 1, -1, 1, 1, 1, 1, -1, -1, 1, - -1, 1, 1, -1, -1, -1, -1, -1, -1, -1, 1, -1, 1, 1, -1, -1, - -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, -1, 1, 1, - 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, - -1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1, 1, - 1, -1, -1, -1, -1, 1, -1, 1, -1, -1, 1, -1, 1, -1, -1, -1, - 1, -1, 1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, 1, 1, 1, - 1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, 1, - 1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, -1, 1, 1, -1, - -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, - -1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, - -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, - 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, -1, 1, -1, - -1, -1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, -1, - -1, -1, 1, -1, -1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, - -1, -1, 1, -1, -1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, - -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, - 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, - 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, -1, 1, -1, -1, - -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, -1, - -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, -1, -1, -1, - 1, -1, 1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, -1, 1, - -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, - -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, 1, 1, -1, -1, - -1, 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, - -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, 1, - 1, 1, -1, 1, -1, 1, -1, 1, 1, -1, -1, 1, 1, -1, -1, -1, - -1, 1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, - -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, - -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, - -1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, 1, -1, -1, - -1, 1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, -1, 1, 1, - 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, -1, - -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, 1, 1, -1, - 1, -1, -1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, 1, - 1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, - 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, 1, -1, -1, -1, - -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, 1, 1, 1, -1, - 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1, 1, -1, - -1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, - -1, -1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, 1, - -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, -1, -1, 1, -1, -1, 1, - 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, - -1, 1, -1, -1, -1, -1, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, - -1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1, - 1, 1, -1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, -1, - -1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, 1, -1, - 1, 1, 1, 1, -1, -1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, - -1, -1, 1, 1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, - 1, 1, 1, -1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -} +// --- GF tables --- -// GF(2^8) tables with primitive polynomial 0x11d. var gfExp = [512]byte{1, 2, 4, 8, 16, 32, 64, 128, 29, 58, 116, 232, 205, 135, 19, 38, 76, 152, 45, 90, 180, 117, 234, 201, 143, 3, 6, 12, 24, 48, 96, 192, 157, 39, 78, 156, 37, 74, 148, 53, 106, 212, 181, 119, 238, 193, 159, 35, 70, 140, 5, 10, 20, 40, 80, 160, 93, 186, 105, 210, 185, 111, 222, 161, 95, 190, 97, 194, 153, 47, 94, 188, 101, 202, 137, 15, 30, 60, 120, 240, 253, 231, 211, 187, 107, 214, 177, 127, 254, 225, 223, 163, 91, 182, 113, 226, 217, 175, 67, 134, 17, 34, 68, 136, 13, 26, 52, 104, 208, 189, 103, 206, 129, 31, 62, 124, 248, 237, 199, 147, 59, 118, 236, 197, 151, 51, 102, 204, 133, 23, 46, 92, 184, 109, 218, 169, 79, 158, 33, 66, 132, 21, 42, 84, 168, 77, 154, 41, 82, 164, 85, 170, 73, 146, 57, 114, 228, 213, 183, 115, 230, 209, 191, 99, 198, 145, 63, 126, 252, 229, 215, 179, 123, 246, 241, 255, 227, 219, 171, 75, 150, 49, 98, 196, 149, 55, 110, 220, 165, 87, 174, 65, 130, 25, 50, 100, 200, 141, 7, 14, 28, 56, 112, 224, 221, 167, 83, 166, 81, 162, 89, 178, 121, 242, 249, 239, 195, 155, 43, 86, 172, 69, 138, 9, 18, 36, 72, 144, 61, 122, 244, 245, 247, 243, 251, 235, 203, 139, 11, 22, 44, 88, 176, 125, 250, 233, 207, 131, 27, 54, 108, 216, 173, 71, 142, 1, 2, 4, 8, 16, 32, 64, 128, 29, 58, 116, 232, 205, 135, 19, 38, 76, 152, 45, 90, 180, 117, 234, 201, 143, 3, 6, 12, 24, 48, 96, 192, 157, 39, 78, 156, 37, 74, 148, 53, 106, 212, 181, 119, 238, 193, 159, 35, 70, 140, 5, 10, 20, 40, 80, 160, 93, 186, 105, 210, 185, 111, 222, 161, 95, 190, 97, 194, 153, 47, 94, 188, 101, 202, 137, 15, 30, 60, 120, 240, 253, 231, 211, 187, 107, 214, 177, 127, 254, 225, 223, 163, 91, 182, 113, 226, 217, 175, 67, 134, 17, 34, 68, 136, 13, 26, 52, 104, 208, 189, 103, 206, 129, 31, 62, 124, 248, 237, 199, 147, 59, 118, 236, 197, 151, 51, 102, 204, 133, 23, 46, 92, 184, 109, 218, 169, 79, 158, 33, 66, 132, 21, 42, 84, 168, 77, 154, 41, 82, 164, 85, 170, 73, 146, 57, 114, 228, 213, 183, 115, 230, 209, 191, 99, 198, 145, 63, 126, 252, 229, 215, 179, 123, 246, 241, 255, 227, 219, 171, 75, 150, 49, 98, 196, 149, 55, 110, 220, 165, 87, 174, 65, 130, 25, 50, 100, 200, 141, 7, 14, 28, 56, 112, 224, 221, 167, 83, 166, 81, 162, 89, 178, 121, 242, 249, 239, 195, 155, 43, 86, 172, 69, 138, 9, 18, 36, 72, 144, 61, 122, 244, 245, 247, 243, 251, 235, 203, 139, 11, 22, 44, 88, 176, 125, 250, 233, 207, 131, 27, 54, 108, 216, 173, 71, 142, 1, 2} var gfLog = [256]byte{0, 0, 1, 25, 2, 50, 26, 198, 3, 223, 51, 238, 27, 104, 199, 75, 4, 100, 224, 14, 52, 141, 239, 129, 28, 193, 105, 248, 200, 8, 76, 113, 5, 138, 101, 47, 225, 36, 15, 33, 53, 147, 142, 218, 240, 18, 130, 69, 29, 181, 194, 125, 106, 39, 249, 185, 201, 154, 9, 120, 77, 228, 114, 166, 6, 191, 139, 98, 102, 221, 48, 253, 226, 152, 37, 179, 16, 145, 34, 136, 54, 208, 148, 206, 143, 150, 219, 189, 241, 210, 19, 92, 131, 56, 70, 64, 30, 66, 182, 163, 195, 72, 126, 110, 107, 58, 40, 84, 250, 133, 186, 61, 202, 94, 155, 159, 10, 21, 121, 43, 78, 212, 229, 172, 115, 243, 167, 87, 7, 112, 192, 247, 140, 128, 99, 13, 103, 74, 222, 237, 49, 197, 254, 24, 227, 165, 153, 119, 38, 184, 180, 124, 17, 68, 146, 217, 35, 32, 137, 46, 55, 63, 209, 91, 149, 188, 207, 205, 144, 135, 151, 178, 220, 252, 190, 97, 242, 86, 211, 171, 20, 42, 93, 158, 132, 60, 57, 83, 71, 109, 65, 162, 31, 45, 67, 216, 183, 123, 164, 118, 196, 23, 73, 236, 127, 12, 111, 246, 108, 161, 59, 82, 41, 157, 85, 170, 251, 96, 134, 177, 187, 204, 62, 90, 203, 89, 95, 176, 156, 169, 160, 81, 11, 245, 22, 235, 122, 117, 44, 215, 79, 174, 213, 233, 230, 231, 173, 232, 116, 214, 244, 234, 168, 80, 88, 175} var rsGen = [9]byte{1, 255, 11, 81, 54, 239, 173, 200, 24} - - -// rsGen is the RS(16,8) generator polynomial coefficients (fcr=0). diff --git a/internal/watermark/watermark_roundtrip_test.go b/internal/watermark/watermark_roundtrip_test.go deleted file mode 100644 index 84fab12..0000000 --- a/internal/watermark/watermark_roundtrip_test.go +++ /dev/null @@ -1,188 +0,0 @@ -package watermark - -import ( - "math" - "sort" - "testing" -) - -// TestRoundTrip verifies the full embed → downsample → phase-search → rotation → RS-decode chain. -func TestRoundTrip(t *testing.T) { - const key = "test-key-42" - const recRate = float64(RecordingRate) // 48000 - const compRate = float64(CompositeRate) // 228000 - const duration = 60.0 // seconds — ~2.7 frames at ChipRate=12kHz - - nRecSamples := int(duration * recRate) - - // === Embed === - emb := NewEmbedder(key) - // No gate — test pure watermark signal - samples := make([]float64, 0, nRecSamples) - - // Drive embedder at CompositeRate, collect at RecordingRate via Bresenham - accum := 0 - var last float64 - for len(samples) < nRecSamples { - last = emb.NextSample() - accum += RecordingRate - if accum >= CompositeRate { - accum -= CompositeRate - samples = append(samples, last) - } - } - - t.Logf("Embedded: %d samples @ %.0f Hz = %.2fs", len(samples), recRate, float64(len(samples))/recRate) - - // RMS check - var rmsAcc float64 - for _, s := range samples { - rmsAcc += s * s - } - rms := math.Sqrt(rmsAcc / float64(len(samples))) - rmsDBFS := 20 * math.Log10(rms+1e-12) - t.Logf("Watermark RMS: %.1f dBFS (expect ~-48)", rmsDBFS) - if rmsDBFS < -52 || rmsDBFS > -44 { - t.Errorf("RMS %.1f dBFS out of expected range [-52, -44]", rmsDBFS) - } - - // === Decode: Phase search === - samplesPerBit := int(float64(PnChips) * recRate / float64(ChipRate)) - t.Logf("samplesPerBit=%d, frameLen=%d", samplesPerBit, samplesPerBit*PayloadBits) - - const coarseStep = 8 - const syncBits = 64 - - bestPhase := 0 - bestMag := 0.0 - for phase := 0; phase < samplesPerBit; phase += coarseStep { - mag := testAvgCorrMag(samples, phase, samplesPerBit, syncBits, recRate) - if mag > bestMag { - bestMag = mag - bestPhase = phase - } - } - fineStart := bestPhase - coarseStep - if fineStart < 0 { fineStart = 0 } - fineEnd := bestPhase + coarseStep - if fineEnd > samplesPerBit { fineEnd = samplesPerBit } - for phase := fineStart; phase < fineEnd; phase++ { - mag := testAvgCorrMag(samples, phase, samplesPerBit, syncBits, recRate) - if mag > bestMag { - bestMag = mag - bestPhase = phase - } - } - t.Logf("Phase search: bestPhase=%d, avgCorr=%.4f", bestPhase, bestMag) - - // Phase should be 0 for clean signal starting at sample 0 - if bestPhase != 0 { - t.Errorf("expected bestPhase=0, got %d", bestPhase) - } - - // === Decode: Extract correlations === - nCompleteBits := (len(samples) - bestPhase) / samplesPerBit - nFrames := nCompleteBits / PayloadBits - if nFrames == 0 { nFrames = 1 } - t.Logf("Complete bits: %d, frames: %d", nCompleteBits, nFrames) - - corrs := make([]float64, PayloadBits) - for i := 0; i < PayloadBits; i++ { - for frame := 0; frame < nFrames; frame++ { - bitGlobal := frame*PayloadBits + i - start := bestPhase + bitGlobal*samplesPerBit - if start+samplesPerBit > len(samples) { break } - corrs[i] += CorrelateAt(samples, start, recRate) - } - } - - // Log correlation stats - var minAbs, maxAbs float64 - for i, c := range corrs { - ac := math.Abs(c) - if i == 0 || ac < minAbs { minAbs = ac } - if ac > maxAbs { maxAbs = ac } - } - t.Logf("Correlation range: min|c|=%.2f, max|c|=%.2f", minAbs, maxAbs) - - // === Decode: Frame sync via rotation === - type decodeResult struct { - rotation int - payload [RsDataBytes]byte - erasures int - } - var best *decodeResult - - for rot := 0; rot < PayloadBits; rot++ { - var recv [RsTotalBytes]byte - confs := make([]float64, PayloadBits) - for i := 0; i < PayloadBits; i++ { - srcBit := (i + rot) % PayloadBits - c := corrs[srcBit] - confs[i] = math.Abs(c) - if c < 0 { - recv[i/8] |= 1 << uint(7-(i%8)) - } - } - - type bitConf struct { idx int; conf float64 } - ranked := make([]bitConf, PayloadBits) - for i := range ranked { ranked[i] = bitConf{i, confs[i]} } - sort.Slice(ranked, func(a, b int) bool { return ranked[a].conf < ranked[b].conf }) - - for nErase := 0; nErase <= RsCheckBytes*8; nErase++ { - erasedBytes := map[int]bool{} - for _, bc := range ranked[:nErase] { - erasedBytes[bc.idx/8] = true - } - if len(erasedBytes) > RsCheckBytes { break } - erasePos := make([]int, 0, len(erasedBytes)) - for pos := range erasedBytes { erasePos = append(erasePos, pos) } - sort.Ints(erasePos) - - payload, ok := RSDecode(recv, erasePos) - if ok { - if best == nil || len(erasePos) < best.erasures { - best = &decodeResult{rotation: rot, payload: payload, erasures: len(erasePos)} - } - break - } - } - if best != nil && best.erasures == 0 { break } - } - - if best == nil { - t.Fatal("RS decode FAILED — no valid rotation found") - } - - t.Logf("Decoded: rotation=%d, erasures=%d, payload=%x", best.rotation, best.erasures, best.payload) - - // Rotation should be 0 for clean signal - if best.rotation != 0 { - t.Errorf("expected rotation=0, got %d", best.rotation) - } - if best.erasures != 0 { - t.Errorf("expected 0 erasures, got %d", best.erasures) - } - - // Key match - if !KeyMatchesPayload(key, best.payload) { - t.Errorf("key %q does NOT match decoded payload %x", key, best.payload) - } else { - t.Logf("Key %q MATCHES", key) - } -} - -func testAvgCorrMag(samples []float64, phase, samplesPerBit, nBits int, recRate float64) float64 { - var total float64 - var count int - for b := 0; b < nBits; b++ { - start := phase + b*samplesPerBit - if start+samplesPerBit > len(samples) { break } - c := CorrelateAt(samples, start, recRate) - total += math.Abs(c) - count++ - } - if count == 0 { return 0 } - return total / float64(count) -}