瀏覽代碼

watermark: replace time-domain PN with STFT-domain spread-spectrum (Kirovski & Malvar 2003)

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
main
Jan 1 月之前
父節點
當前提交
e8c5c2729b
共有 7 個文件被更改,包括 222 次插入765 次删除
  1. +0
    -3
      cmd/fmrtx/main.go
  2. +122
    -130
      cmd/wmdecode/main.go
  3. +9
    -6
      internal/offline/generator.go
  4. +4
    -2
      internal/watermark/stft_roundtrip_test.go
  5. +20
    -16
      internal/watermark/stft_watermark.go
  6. +67
    -420
      internal/watermark/watermark.go
  7. +0
    -188
      internal/watermark/watermark_roundtrip_test.go

+ 0
- 3
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


+ 122
- 130
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 <file.wav> [key ...]
// wmdecode <file.wav> Extract fingerprint
// wmdecode <file.wav> [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 {


+ 9
- 6
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
}



+ 4
- 2
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))


+ 20
- 16
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)


+ 67
- 420
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).

+ 0
- 188
internal/watermark/watermark_roundtrip_test.go 查看文件

@@ -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)
}

Loading…
取消
儲存