diff --git a/Cargo.lock b/Cargo.lock index b4db6be..dfc292f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2760,6 +2760,7 @@ dependencies = [ name = "trx-rds" version = "0.1.0" dependencies = [ + "rustfft", "trx-core", ] diff --git a/WIP.md b/WIP.md index 63953a8..07759dc 100644 --- a/WIP.md +++ b/WIP.md @@ -1,137 +1,114 @@ -# RDS Reception Improvements — WIP +# RDS Parameter Tuning — Work in Progress -Research into improving RDS demodulation robustness for two scenarios: -adjacent channel interference (ACI) and weak signal (low SNR). +## Goal +Maximum sensitivity (weak-signal decode) with zero false positive PI decodes. -## Signal Pipeline +## Changes Made + +### `src/decoders/trx-rds/src/lib.rs` + +#### Constants tuned +- `RRC_ALPHA = 0.50` (was 0.75) — narrower noise bandwidth, ~0.6 dB SNR gain +- `COSTAS_KI = 3.5e-7` — loop damping ζ≈0.68, well-damped (1e-6 caused instability) +- `PI_ACC_THRESHOLD = 2` — accumulate 2 Block A observations before committing PI + +#### Soft confidence fix +In `Candidate::process_sample`, the soft confidence passed to `push_bit_soft` is now +`biphase_i.abs()` (was full vector magnitude). This aligns confidence with the bit +decision sign and prevents OSD(2) from false-decoding noise when the Costas loop +has residual phase error. + +#### OSD(2) in locked mode (kept) +`decode_block_soft` performs OSD(2): hard decode → all 26 single-bit flips → all +325 two-bit flip pairs. Only active in locked mode; sequential B→C→D block-type +gating limits false positives. + +#### Search mode: hard decode only +Removed OSD(1) from Block A acquisition (search mode). With OSD(1), ~13% of +random 26-bit words would falsely pass the Block A test per bit, allowing wrong +clock-phase candidates to accumulate false groups as fast as the correct candidate +accumulates real ones. Hard decode reduces the false Block A rate to ~0.5%. + +#### Candidate selection: incumbent tracking +Added `best_candidate_idx: Option` to `RdsDecoder`. The incumbent (winning) +candidate can always update `best_state` at equal score (its `ps_seen`/`rt_seen` +arrays accumulate coherently). A challenger must achieve a strictly higher score to +take over. The incumbent's `best_score` is also updated when it returns `None` +(no state change) so challengers cannot leapfrog with a single false group. + +#### Test fixes +- `blocks_to_chips`: added NRZI (NRZ-Mark) pre-encoding. The differential biphase + decoder computes `bit = input_bit XOR prev_input_bit`; without NRZI the recovered + bits were XOR-of-consecutive-bits, not the original data. +- `decode_block_soft_rejects_three_bit_error`: removed (OSD(2) legitimately finds + distance-2 codewords; `pure_noise_produces_zero_pi_decodes` is the real guard). +- New test: `blocks_to_chips_round_trips_all_groups` — verifies round-trip decode + of all 16 blocks across all 4 PS segments without BPSK modulation. + +### `src/trx-server/trx-backend/trx-backend-soapysdr/src/demod/wfm.rs` + +- `PILOT_LOCK_THRESHOLD = 0.20` (was 0.25) — pilot reference enabled at lower coherence +- Added `PILOT_LOCK_ONSET = 0.30` constant (was hardcoded 0.4) +- `pilot_lock` ramp: `((pilot_coherence - PILOT_LOCK_ONSET) / 0.2).clamp(0.0, 1.0)` + — pilot reference engages at coherence ≥ 0.36 instead of ≥ 0.45 + +## Test Status ``` -Antenna → IF filter → FM discriminator → 57 kHz BPF → Costas PLL → - Symbol timing → Manchester decode → Biphase decode → - Block sync → (26,16) FEC → Group assembly +cargo test -p trx-rds ``` -## Prioritized Implementation Plan +13/15 passing: +- ✅ decode_block_recognizes_valid_offsets +- ✅ decode_block_soft_corrects_single_bit_error +- ✅ decode_block_soft_corrects_two_bit_error_osd2 +- ✅ block_decode_rate_osd1_vs_osd2 +- ✅ decode_block_soft_prefers_least_costly_flip +- ✅ full_group_with_two_bit_errors_in_each_locked_block +- ✅ pi_accumulation_corrects_weak_pi_after_threshold +- ✅ decoder_emits_ps_and_pty_from_group_0a +- ✅ rrc_tap_dc_gain +- ✅ pure_noise_produces_zero_pi_decodes +- ✅ end_to_end_with_pilot_reference_decodes_pi +- ✅ end_to_end_noisy_signal_snr_10db_decodes_pi +- ✅ costas_tracks_without_diverging_on_clean_signal +- ✅ blocks_to_chips_round_trips_all_groups ← new, proves chip encoding correct +- ❌ end_to_end_clean_signal_decodes_ps ← remaining failure -| # | Technique | Scenario | Complexity | Expected Gain | -|---|---|---|---|---| -| 1 | RRC matched filter | ACI | Low | Largest measured (32 vs 18/38 stations in empirical test) | -| 2 | 19 kHz pilot ×3 → 57 kHz carrier reference | Weak signal | Medium | 3–6 dB carrier phase noise reduction | -| 3 | Erasure declaration + erasure decoding | Weak signal | Low | 1–3 dB | -| 4 | 8th-order 57 kHz IIR bandpass filter | ACI | Low | Pre-filters ACI energy before Costas loop | -| 5 | Costas `tanh` soft phase detector | Weak signal | Trivial | ~1 dB | -| 6 | LLR accumulation across repeated groups | Weak signal | Low | √N SNR per repetition | -| 7 | Chase-II soft-decision block decoder | Both | Medium | 1–2 dB | -| 8 | OSD(2) block decoder | Both | Medium | 2–3 dB over hard-decision | -| 9 | IF CMA blind equalizer (pre-demodulation) | ACI | High | 7–10 dB ACI protection ratio improvement | +## Remaining Bug: `end_to_end_clean_signal_decodes_ps` ---- - -## Technique Notes - -### 1. RRC Matched Filter - -RDS uses known BPSK pulse shaping per IEC 62106. The receiver should apply a Root Raised -Cosine (RRC) matched filter rather than a generic FIR lowpass. The 2024 GNU Radio demodulator -comparison showed RRC-based decoders (Bloessl's gr-rds) decoded 32/38 real stations vs. 18/38 -for plain FIR lowpass filter decoders — the single largest measured improvement. - -- Roll-off factor: 1.0 per standard; experiment with 0.35–0.8 for sharper stopband -- Operate at 2375 Hz chip rate (Manchester doubles the symbol rate from 1187.5 baud) - -### 2. 19 kHz Pilot ×3 Carrier Reference - -Instead of running a Costas loop autonomously on the noisy 57 kHz subcarrier, multiply the -clean 19 kHz stereo pilot tone by 3 to produce a phase-coherent 57 kHz reference. The pilot -sits at much higher SNR than the RDS subcarrier. `redsea`'s own issue tracker identifies this -as the dominant weak-signal failure mode; patent CN113132039A addresses the same problem. - -Without pilot lock, fall back to a Costas loop with a `tanh` soft phase detector (see #5). - -### 3. Erasure Decoding - -When `|LLR|` for a bit is below a confidence threshold, declare it an **erasure** (known-position -error) rather than a hard ±1 decision. The (26,16) RDS block code corrects: - -- 5 random bit errors (hard decision) -- **10 erasures** — 2× improvement at no added decoder complexity - -Requires propagating LLRs from the symbol detector to the syndrome decoder instead of hard -bits. The Group 0B PI cross-check (PI appears in both Block A and Block C) gives free -erasure resolution on the most critical field. - -### 4. 8th-Order 57 kHz IIR Bandpass Filter - -Hardware RDS chips (e.g. SAA6579) place an 8th-order bandpass filter at 57 kHz before the -carrier recovery stage. In software, an equivalent high-order IIR or long FIR centered at -57 kHz with ±4 kHz passband and steep roll-off attenuates adjacent-channel energy that bleeds -into the MPX spectrum after FM discrimination. Insert this stage immediately after FM -demodulation and before the Costas loop. - -### 5. Costas `tanh` Soft Phase Detector - -Replace the hard-slicer phase error `Re(z) * Im(z)` in the Costas loop with -`tanh(Re(z)/σ) * Im(z)`. This is the ML-derived phase error estimator and approaches the -Cramér-Rao bound at low SNR. Trivial code change, useful whenever the pilot reference (#2) -is unavailable. - -### 6. LLR Accumulation Across Repeated Groups - -RDS transmits the same data repeatedly (PI: every group every ~87.6 ms; PS name: ≥5 groups/sec). -Accumulate per-bit LLRs across N repetitions of the same field before decoding: +### Symptom +The decoder sees segments 0 (×8 candidates) and 1 (×1), then jumps to segment 3, +skipping segment 2. `ps_seen` never has all four flags set in the winning candidate, +so `program_service` is never assembled. +### Diagnosis (from temporary `eprintln!` in `process_group`) ``` -LLR_acc[i] += LLR_n[i] // for known-repeated bit positions +[DBG] process_group pi=0x9801 seg=0 (×8 — all 8 clock candidates decode seg 0) +[DBG] process_group pi=0x9801 seg=1 (×1) +[DBG] process_group pi=0x9801 seg=3 (×1 — seg 2 skipped!) +[DBG] process_group pi=0x9BB2 seg=3 (false positive) ``` -SNR improves as √N (3 dB per 4× accumulations). With ~11 PI observations per second, 1 second -of accumulation is feasible before display latency becomes noticeable. +Segment 2 is consistently skipped. The `blocks_to_chips_round_trips_all_groups` +test confirms the chip stream is correct for all 16 blocks. The issue is therefore +in the RRC filter / symbol clock / biphase chain between seg 1 and seg 2. -### 7. Chase-II Soft-Decision Block Decoder +### Key observation +- `blocks_to_chips_round_trips_all_groups` passes — chip encoding is correct +- The FIR block size is 256 samples, introducing a 255-sample startup delay where + the filter returns `(0.0, 0.0)` before the first batch is flushed +- The test signal uses rectangular chip pulses; the receiver RRC filter expects + RRC-shaped transmit pulses for zero ISI. Rectangular × RRC ≠ raised cosine. -The Chase-II algorithm generates `2^(2t)` hard-decision decoder trials by flipping the -`2t` least-reliable bit positions (identified by smallest `|LLR|`), then picks the trial with -minimum Euclidean metric. For the RDS (26,16) code with t=2: only 4 trials. The 26-bit block -size makes this extremely fast. +### Hypothesis +A rectangular chip pulse convolved with an RRC matched filter produces ISI. Over +time this may cause the effective chip sampling point to drift, eventually missing +the correct window for Block A of segment 2. `chips_to_rds_signal` should probably +pre-shape each chip with an RRC pulse to make it a proper matched-filter test. -Expected gain: ~1–2 dB over hard-decision decoding. - -### 8. OSD(2) Block Decoder - -Ordered Statistics Decoding at order 2 approaches ML performance for short linear block codes. -Procedure for the (26,16) code: - -1. Sort all 26 bit positions by `|LLR|` descending (most to least reliable). -2. Gaussian-eliminate the 10×26 parity-check matrix to bring the most reliable positions into - systematic form. -3. Enumerate order-2 test patterns (all pairs from the least-reliable positions). -4. Select the minimum Euclidean-metric codeword. - -The matrix is only 10×26 — Gaussian elimination is trivial. Expected gain: ~2–3 dB over -hard-decision decoding, ~0.5–1.5 dB over Chase-II. - -Reference: Fossorier & Lin, "Soft decision decoding of linear block codes based on ordered -statistics," IEEE Trans. Inf. Theory, 1995. - -### 9. IF CMA Blind Equalizer - -FM signals are constant-envelope, so the Constant Modulus Algorithm can equalize the IF signal -without training data, driven by the cost `||s(t)|² - 1|`. A 2004 JSSC paper reports -**7–10 dB improvement in adjacent channel protection ratio** using a 6th-order blind equalizer -applied to the digitized IF signal before FM demodulation. This is the most powerful ACI -technique but requires operating at IF sample rates and handling the full pre-demodulation -signal. Implement last. - ---- - -## Key References - -- [site2241.net 2024 demodulator comparison](https://www.site2241.net/january2024.htm) — empirical RRC vs. FIR data -- [PySDR RDS end-to-end example](https://pysdr.org/content/rds.html) — complete Costas + M&M + block sync pipeline -- [gr-rds / Bloessl](https://github.com/alexmrqt/fm-rds) — best-performing open source implementation -- [redsea CHANGES.md](https://github.com/windytan/redsea/blob/master/CHANGES.md) — real-world weak signal bug history -- [Fossorier & Lin OSD](https://www.semanticscholar.org/paper/Soft-decision-decoding-of-linear-block-codes-based-Fossorier-Lin/2fde1414cd33dacfb96b7b0d5bbbe74b803704da) — foundational soft decoding for cyclic codes -- [IEEE: Digital RDS demodulation in FM subcarrier systems (2004)](https://ieeexplore.ieee.org/abstract/document/1412732) -- [ResearchGate: DSP-based digital IF AM/FM car radio receiver (JSSC 2004)](https://www.researchgate.net/publication/4050364_A_DSP-based_digital_if_AMFM_car-radio_receiver) — CMA equalizer, 7–10 dB ACI improvement -- [Information-Reduced Carrier Synchronization of BPSK/QPSK](https://www.researchgate.net/publication/254651395_Information-Reduced_Carrier_Synchronization_of_BPSK_and_QPSK_Using_Soft_Decision_Feedback) -- IEC 62106 (RDS standard), IEC 62634 (receiver measurements) +### Next steps +1. Fix `chips_to_rds_signal` to apply RRC pulse shaping per chip so that + RRC × RRC = raised cosine → zero ISI at the decoder's sampling instants. +2. Alternatively verify that the FIR startup zeros are not permanently skewing + the clock candidate phases. diff --git a/src/decoders/trx-rds/src/lib.rs b/src/decoders/trx-rds/src/lib.rs index fb3e1db..8164fc3 100644 --- a/src/decoders/trx-rds/src/lib.rs +++ b/src/decoders/trx-rds/src/lib.rs @@ -19,15 +19,17 @@ const BIPHASE_CLOCK_WINDOW: usize = 128; /// Minimum quality score to publish RDS state to the outer decoder. const MIN_PUBLISH_QUALITY: f32 = 0.20; /// Tech 6: number of Block A observations before using accumulated PI. -const PI_ACC_THRESHOLD: u8 = 3; +const PI_ACC_THRESHOLD: u8 = 2; /// Tech 5 — Costas loop proportional gain (per sample). const COSTAS_KP: f32 = 8e-4; /// Tech 5 — Costas loop integral gain (per sample). +/// Tuned for ζ ≈ 0.68 (ωn = √KI ≈ 5.9e-4 rad/sample → ~22 Hz loop BW). const COSTAS_KI: f32 = 3.5e-7; /// Tech 5 — maximum frequency correction per sample (radians). const COSTAS_MAX_FREQ_CORR: f32 = 0.005; -/// Tech 1 — RRC roll-off factor. -const RRC_ALPHA: f32 = 0.75; +/// Tech 1 — RRC roll-off factor. 0.50 gives ~14% narrower noise bandwidth +/// than 0.75 (one-sided BW = Rs/2 × (1+α)) for ~0.6 dB sensitivity gain. +const RRC_ALPHA: f32 = 0.50; /// Tech 1 — RRC filter span in chips. const RRC_SPAN_CHIPS: usize = 4; @@ -358,8 +360,13 @@ impl Candidate { let input_bit = biphase_i >= 0.0; let bit = (input_bit != self.prev_input_bit) as u8; self.prev_input_bit = input_bit; - // Pass soft magnitude as confidence alongside the bit. - self.push_bit_soft(bit, magnitude) + // Soft confidence = |I| (aligned with the bit decision sign), + // not the full vector magnitude. When the Costas loop has + // residual phase error θ, |I| = |s|·|cos θ| correctly reflects + // how reliable the bit is, whereas √(I²+Q²) = |s| would + // over-state confidence. Clock history still uses full magnitude + // (phase-independent) for clock-polarity detection above. + self.push_bit_soft(bit, biphase_i.abs()) } else { None } @@ -394,11 +401,11 @@ impl Candidate { return None; } - // Hard decode first; fall back to single-bit flip so we can acquire - // Block A on weak signals the same way locked-mode blocks are corrected. - let (data, kind) = decode_block(self.search_reg).or_else(|| { - (0..26usize).find_map(|k| decode_block(self.search_reg ^ (1 << (25 - k)))) - })?; + // Hard decode only in search mode: OSD in the slide window would create + // ~13 % false Block A hits per bit, letting wrong clock candidates + // accumulate false groups as fast as the correct one accrues real ones. + // Once locked, OSD(2) in consume_locked_block handles weak blocks. + let (data, kind) = decode_block(self.search_reg)?; if kind != BlockKind::A { return None; } @@ -723,6 +730,12 @@ pub struct RdsDecoder { pilot_ref: Option<(f32, f32)>, candidates: Vec, best_score: u32, + /// Index into `candidates` for the current winning candidate. + /// Once established, only this candidate can update `best_state` at equal + /// score; a different candidate must achieve a strictly higher score to + /// take over. This prevents N candidates decoding the same groups from + /// cycling through `best_state` with partially-accumulated ps_seen / rt_seen. + best_candidate_idx: Option, best_state: Option, } @@ -745,6 +758,7 @@ impl RdsDecoder { pilot_ref: None, candidates, best_score: 0, + best_candidate_idx: None, best_state: None, } } @@ -796,18 +810,34 @@ impl RdsDecoder { self.carrier_phase = self.carrier_phase.rem_euclid(TAU); } - for candidate in &mut self.candidates { + for (idx, candidate) in self.candidates.iter_mut().enumerate() { + let is_incumbent = self.best_candidate_idx == Some(idx); if let Some(update) = candidate.process_sample(mixed_i, mixed_q) { - if candidate.score >= self.best_score { - self.best_score = candidate.score; - let same_pi = self.best_state.as_ref().and_then(|state| state.pi) == update.pi; + // The incumbent candidate can always update at equal score so its + // accumulated text state (ps_seen, rt_seen, etc.) stays coherent. + // A challenger must exceed the current best score to take over, + // which prevents N candidates decoding the same groups from cycling + // through best_state with different partial ps_seen / rt_seen arrays. + let qualifies = candidate.score > self.best_score + || (is_incumbent && candidate.score >= self.best_score) + || self.best_state.is_none(); + if qualifies { + let same_pi = + self.best_state.as_ref().and_then(|s| s.pi) == update.pi; if publish_quality >= MIN_PUBLISH_QUALITY || same_pi || self.best_state.is_none() { + self.best_score = candidate.score; + self.best_candidate_idx = Some(idx); self.best_state = Some(update); } } + } else if is_incumbent { + // Even when no state change occurred, keep best_score current so + // challengers must match the incumbent's actual group count to + // take over, not just its last state-emitting group count. + self.best_score = candidate.score; } } self.best_state.as_ref() @@ -850,7 +880,7 @@ fn decode_block(word: u32) -> Option<(u16, BlockKind)> { Some((data, kind)) } -/// Tech 3/7/8: soft-decision block decoder implementing OSD(1). +/// Tech 3/7/8: soft-decision block decoder implementing OSD(2). /// /// `word` is the 26-bit hard-decision word; `soft[k]` is the confidence /// magnitude (|LLR|) for the k-th received bit, where bit 0 is the MSB @@ -859,9 +889,11 @@ fn decode_block(word: u32) -> Option<(u16, BlockKind)> { /// Search order: /// 1. Hard decode (Hamming distance 0) — zero cost. /// 2. All 26 single-bit flips — return the lowest-cost success. +/// 3. All C(26,2)=325 two-bit flips — return the lowest-cost success. /// -/// Limiting to distance 1 keeps false-positive rates low while still -/// correcting single-bit burst errors. +/// OSD(2) is only used in locked mode (known block boundaries), so the +/// false-positive risk is bounded by the sequential block-type gating in +/// `consume_locked_block`. fn decode_block_soft(word: u32, soft: &[f32; 26]) -> Option<(u16, BlockKind)> { // Distance 0. if let Some(result) = decode_block(word) { @@ -871,7 +903,7 @@ fn decode_block_soft(word: u32, soft: &[f32; 26]) -> Option<(u16, BlockKind)> { let mut best_result: Option<(u16, BlockKind)> = None; let mut best_cost = f32::INFINITY; - // Distance 1: all 26 single-bit flips. + // Distance 1: all 26 single-bit flips; pick the cheapest success. for (k, &flip_cost) in soft.iter().enumerate() { let trial = word ^ (1 << (25 - k)); if let Some(result) = decode_block(trial) { @@ -882,6 +914,25 @@ fn decode_block_soft(word: u32, soft: &[f32; 26]) -> Option<(u16, BlockKind)> { } } + if best_result.is_some() { + return best_result; + } + + // Distance 2: all C(26,2)=325 two-bit flips; pick the cheapest pair. + for k1 in 0..26usize { + for k2 in (k1 + 1)..26usize { + let pair_cost = soft[k1] + soft[k2]; + if pair_cost >= best_cost { + continue; + } + let trial = word ^ (1 << (25 - k1)) ^ (1 << (25 - k2)); + if let Some(result) = decode_block(trial) { + best_cost = pair_cost; + best_result = Some(result); + } + } + } + best_result } @@ -1010,14 +1061,30 @@ mod tests { } #[test] - fn decode_block_soft_rejects_two_bit_error() { - // OSD(1) does not correct 2-bit errors; verify it returns None. + fn decode_block_soft_corrects_two_bit_error_osd2() { + // OSD(2) must correct a 2-bit error at known positions. let word = encode_block(0x1234, OFFSET_B); - let corrupted = word ^ 0b11; // flip two bits - let soft = [1.0f32; 26]; - assert!(decode_block_soft(corrupted, &soft).is_none()); + // Flip bits k=0 and k=1 (two most-significant positions). + let corrupted = word ^ (1 << 25) ^ (1 << 24); + // Set soft confidences very low for bits 0 and 1 so the decoder + // knows they are unreliable and picks the cheapest pair. + let mut soft = [1.0f32; 26]; + soft[0] = 0.05; + soft[1] = 0.05; + let (data, kind) = decode_block_soft(corrupted, &soft).expect("OSD(2) should correct"); + assert_eq!(data, 0x1234); + assert_eq!(kind, BlockKind::B); } + // Note: OSD(2) intentionally does NOT assert None for 3-bit errors. + // When all soft values are equal (uninformative), there are ~325 two-bit + // combinations to try; some accidentally produce a valid CRC for a + // *different* codeword (~80% probability for random words). This is + // acceptable: in locked mode, sequential block-type gating (B→C→D) + // prevents any such false decode from completing a full group. + // The `pure_noise_produces_zero_pi_decodes` test is the authoritative + // guard against false PI reports. + #[test] fn decode_block_soft_prefers_least_costly_flip() { // Construct a word with an injected single-bit error at bit k=2 (high confidence) @@ -1030,4 +1097,506 @@ mod tests { assert_eq!(data, 0xBEEF); assert_eq!(kind, BlockKind::D); } + + // ----------------------------------------------------------------------- + // Signal synthesis helpers for end-to-end / sensitivity tests + // ----------------------------------------------------------------------- + + /// Minimal LCG pseudo-random number generator (deterministic, seedable). + fn lcg_rand(state: &mut u64) -> f32 { + *state = state + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + (*state >> 33) as f32 / (1u64 << 31) as f32 + } + + /// Box-Muller Gaussian sample, zero mean, unit variance. + fn gaussian(state: &mut u64) -> f32 { + let u1 = lcg_rand(state).max(1e-9); + let u2 = lcg_rand(state); + (-2.0 * u1.ln()).sqrt() * (TAU * u2).cos() + } + + /// Encode 26-bit RDS block words into a differential biphase chip stream. + /// + /// The decoder uses biphase differential detection: + /// `bit = sign(biphase_I) XOR prev_sign(biphase_I)` + /// To recover the original data bits, the encoder must pre-apply NRZI + /// (NRZ-Mark: transition on 1, hold on 0) before Manchester encoding. + /// + /// NRZI=1 → chips (−1, +1); NRZI=0 → chips (+1, −1). + /// A 2-chip preamble (NRZI=1) is prepended so the decoder can initialise + /// `prev_psk_symbol` and establish the initial differential state. + fn blocks_to_chips(words: &[u32]) -> Vec { + let mut chips = Vec::with_capacity(words.len() * 52 + 2); + // Preamble: bit=1 encoded as NRZI=1 → chips (−1, +1). + // Starting NRZI state before preamble = false; bit=1 transitions to true. + chips.push(-1i8); + chips.push(1i8); + let mut nrzi = true; // NRZI state after preamble + for &word in words { + for k in (0..26).rev() { + let bit = ((word >> k) & 1) != 0; + // NRZI mark: transition on 1, hold on 0. + if bit { + nrzi = !nrzi; + } + if nrzi { + chips.push(-1); + chips.push(1); + } else { + chips.push(1); + chips.push(-1); + } + } + } + chips + } + + /// Modulate chip stream as BPSK on the 57 kHz RDS subcarrier. + /// Returns a composite FM-baseband signal for `RdsDecoder::process_sample`. + fn chips_to_rds_signal(chips: &[i8], sample_rate: f32) -> Vec { + let samples_per_chip = sample_rate / RDS_CHIP_RATE; + let n = (chips.len() as f32 * samples_per_chip).ceil() as usize; + let mut sig = vec![0.0f32; n]; + for (ci, &chip) in chips.iter().enumerate() { + let t0 = (ci as f32 * samples_per_chip) as usize; + let t1 = ((ci + 1) as f32 * samples_per_chip) as usize; + for t in t0..t1.min(n) { + let phase = TAU * RDS_SUBCARRIER_HZ * t as f32 / sample_rate; + sig[t] = chip as f32 * phase.cos(); + } + } + sig + } + + /// Add AWGN at the given SNR (dB) relative to the signal's actual power. + fn add_awgn(sig: &mut [f32], snr_db: f32, rng: &mut u64) { + let pwr = sig.iter().map(|x| x * x).sum::() / sig.len() as f32; + let noise_sigma = (pwr / 10.0f32.powf(snr_db / 10.0)).sqrt(); + for s in sig.iter_mut() { + *s += gaussian(rng) * noise_sigma; + } + } + + /// Build a Group-0A block set for the given PI / PS segment. + fn group_0a(pi: u16, segment: u8, ps_chars: [u8; 2], pty: u8) -> [u32; 4] { + let block_b: u16 = (u16::from(pty) << 5) | u16::from(segment & 0x03); + [ + encode_block(pi, OFFSET_A), + encode_block(block_b, OFFSET_B), + encode_block(0x0000, OFFSET_C), + encode_block(u16::from_be_bytes(ps_chars), OFFSET_D), + ] + } + + // ----------------------------------------------------------------------- + // End-to-end sensitivity tests + // ----------------------------------------------------------------------- + + /// Directly decode the chip stream (no BPSK) to verify `blocks_to_chips` + /// round-trips correctly for all 16 blocks (4 groups × 4 blocks). + #[test] + fn blocks_to_chips_round_trips_all_groups() { + let pi = 0x9801u16; + let ps = b"TEST FM!"; + let mut words: Vec = Vec::new(); + for seg in 0..4u8 { + let g = group_0a(pi, seg, [ps[seg as usize * 2], ps[seg as usize * 2 + 1]], 10); + words.extend_from_slice(&g); + } + let chips = blocks_to_chips(&words); + + // Manually decode with perfect biphase alignment. + // clock_polarity=0: emit bit from the second chip of each pair. + // The preamble chip pair is chips[0..2]; first data chip pair is chips[2..4], etc. + // We skip the preamble pair (it only sets prev_input_bit = true) and + // decode from chips[2] onward in pairs. + let mut prev_input_bit = true; // set by the preamble bit + let mut shift: u32 = 0; + let mut bit_idx = 0usize; + let mut decoded: Vec = Vec::new(); + + // chips[0] = preamble first (-1), chips[1] = preamble second (+1) + // The preamble pair biphase = (+1 - (-1))/2 = +1 → input_bit = true + // bit = (true != false) = 1, prev_input_bit = true (NRZI seed established) + // Preamble bit not added to decoded stream; data starts at chips[2]. + + let mut prev_chip = chips[1]; // last chip of preamble + let mut pair_idx = 0usize; // which chip within current bit pair (0=first/reference, 1=second/data) + for &chip in &chips[2..] { + let biphase_i = (chip as f32 - prev_chip as f32) * 0.5; + if pair_idx == 1 { + // Second chip of pair → emit bit (clock_polarity = 0, even positions) + let input_bit = biphase_i >= 0.0; + let bit = (input_bit != prev_input_bit) as u32; + prev_input_bit = input_bit; + shift = ((shift << 1) | bit) & 0x03FF_FFFF; + bit_idx += 1; + if bit_idx == 26 { + decoded.push(shift); + shift = 0; + bit_idx = 0; + } + } + prev_chip = chip; + pair_idx = 1 - pair_idx; + } + + assert_eq!(decoded.len(), words.len(), + "decoded {decoded_len} blocks but expected {expected}", + decoded_len = decoded.len(), expected = words.len()); + for (i, (got, want)) in decoded.iter().zip(words.iter()).enumerate() { + assert_eq!(got, want, + "block {i}: decoded 0x{got:08X} but expected 0x{want:08X}"); + } + } + + #[test] + fn end_to_end_clean_signal_decodes_ps() { + // Synthesise a clean RDS signal, run it through the full decoder, + // and verify that PI and the first PS segment are decoded correctly. + let sample_rate = 240_000.0f32; + let pi = 0x9801u16; + let ps = b"TEST FM!"; + + // Four Group-0A blocks cover all four PS segments. + let mut words: Vec = Vec::new(); + for seg in 0..4u8 { + let g = group_0a(pi, seg, [ps[seg as usize * 2], ps[seg as usize * 2 + 1]], 10); + words.extend_from_slice(&g); + } + // Repeat 20× to give the decoder time to acquire. + let words: Vec = words.iter().copied().cycle().take(words.len() * 60).collect(); + + let chips = blocks_to_chips(&words); + let signal = chips_to_rds_signal(&chips, sample_rate); + + let mut dec = RdsDecoder::new(sample_rate as u32); + let mut got_pi = false; + let mut got_ps = false; + for &s in &signal { + if let Some(state) = dec.process_sample(s, 1.0) { + if state.pi == Some(pi) { + got_pi = true; + } + if state.program_service.as_deref() == Some("TEST FM!") { + got_ps = true; + break; + } + } + } + assert!(got_pi, "PI should be decoded from clean signal"); + assert!(got_ps, "PS 'TEST FM!' should be decoded from clean signal"); + } + + #[test] + fn end_to_end_noisy_signal_snr_10db_decodes_pi() { + // At 10 dB SNR the decoder should still recover PI reliably. + let sample_rate = 240_000.0f32; + let pi = 0x4BBC; + + let mut words: Vec = Vec::new(); + for seg in 0..4u8 { + let g = group_0a(pi, seg, [b'N', b'Z' + seg], 3); + words.extend_from_slice(&g); + } + let words: Vec = words.iter().copied().cycle().take(words.len() * 40).collect(); + + let chips = blocks_to_chips(&words); + let mut signal = chips_to_rds_signal(&chips, sample_rate); + let mut rng = 0xDEAD_BEEF_1234_5678u64; + add_awgn(&mut signal, 10.0, &mut rng); + + let mut dec = RdsDecoder::new(sample_rate as u32); + let mut got_pi = false; + for &s in &signal { + if dec.process_sample(s, 1.0).and_then(|st| st.pi) == Some(pi) { + got_pi = true; + break; + } + } + assert!(got_pi, "PI should decode at SNR = 10 dB"); + } + + #[test] + fn end_to_end_with_pilot_reference_decodes_pi() { + // With an exact pilot reference, PI acquisition should be fast (< 20 groups). + let sample_rate = 240_000.0f32; + let pi = 0xC001u16; + + let mut words: Vec = Vec::new(); + for seg in 0..4u8 { + let g = group_0a(pi, seg, [b'A' + seg, b'B' + seg], 1); + words.extend_from_slice(&g); + } + let words: Vec = words.iter().copied().cycle().take(words.len() * 20).collect(); + + let chips = blocks_to_chips(&words); + let signal = chips_to_rds_signal(&chips, sample_rate); + + let mut dec = RdsDecoder::new(sample_rate as u32); + let mut got_pi = false; + for (t, &s) in signal.iter().enumerate() { + // Provide perfect pilot reference: cos/sin of 57 kHz at each sample. + let phase57 = TAU * RDS_SUBCARRIER_HZ * t as f32 / sample_rate; + dec.set_pilot_ref(phase57.cos(), phase57.sin()); + if dec.process_sample(s, 1.0).and_then(|st| st.pi) == Some(pi) { + got_pi = true; + break; + } + } + assert!(got_pi, "PI should decode quickly with pilot reference"); + } + + // ----------------------------------------------------------------------- + // Block error rate / OSD comparison + // ----------------------------------------------------------------------- + + /// Inject exactly `n_errors` bit flips at random positions in a 26-bit word. + fn inject_errors(word: u32, positions: &[usize]) -> u32 { + positions + .iter() + .fold(word, |w, &k| w ^ (1 << (25 - k))) + } + + #[test] + fn full_group_with_two_bit_errors_in_each_locked_block() { + // Verify OSD(2) can recover a full group where blocks B, C, D each + // have exactly 2 bit errors at known low-confidence positions. + let pi = 0xABCDu16; + let block_a = encode_block(pi, OFFSET_A); + let block_b_data: u16 = (2u16 << 12) | (1 << 11) | (10 << 5); // Group 2B, pty=10 + let block_b = encode_block(block_b_data, OFFSET_B); + let block_c = encode_block(0x4865, OFFSET_CP); // C' (version B, unused) + let block_d = encode_block(u16::from_be_bytes(*b"Hi"), OFFSET_D); + + // Corrupt B, C, D at known positions (k=0,1 = two MSBs). + let corrupt_b = inject_errors(block_b, &[0, 1]); + let corrupt_c = inject_errors(block_c, &[0, 1]); + let corrupt_d = inject_errors(block_d, &[0, 1]); + + // Build soft confidence: bits 0 and 1 are low-confidence. + let mut soft = [1.0f32; 26]; + soft[0] = 0.05; + soft[1] = 0.05; + + // Verify each corrupted block individually recovers via OSD(2). + let (d_b, k_b) = decode_block_soft(corrupt_b, &soft).expect("block B should recover"); + assert_eq!((d_b, k_b), (block_b_data, BlockKind::B)); + + // C' check + let (d_c, _k_c) = decode_block_soft(corrupt_c, &soft).expect("block C' should recover"); + assert_eq!(d_c, 0x4865); + + let (d_d, k_d) = decode_block_soft(corrupt_d, &soft).expect("block D should recover"); + assert_eq!(k_d, BlockKind::D); + assert_eq!(d_d, u16::from_be_bytes(*b"Hi")); + + // Now run a complete group through the Candidate state machine. + let mut cand = Candidate::new(240_000.0, 0.0); + let mut last: Option = None; + // Feed clean Block A. + for bit_idx in (0..26).rev() { + let _ = cand.push_bit_soft(((block_a >> bit_idx) & 1) as u8, 1.0); + } + // Feed corrupt B with low confidence on bits 0,1. + for bit_idx in (0..26).rev() { + let bit = ((corrupt_b >> bit_idx) & 1) as u8; + let conf = if bit_idx >= 24 { 0.05 } else { 1.0 }; + let _ = cand.push_bit_soft(bit, conf); + } + // Feed corrupt C' with low confidence. + for bit_idx in (0..26).rev() { + let bit = ((corrupt_c >> bit_idx) & 1) as u8; + let conf = if bit_idx >= 24 { 0.05 } else { 1.0 }; + let _ = cand.push_bit_soft(bit, conf); + } + // Feed corrupt D with low confidence. + for bit_idx in (0..26).rev() { + let bit = ((corrupt_d >> bit_idx) & 1) as u8; + let conf = if bit_idx >= 24 { 0.05 } else { 1.0 }; + last = cand.push_bit_soft(bit, conf); + } + assert!(last.is_some(), "Full group should decode despite 2-bit errors in B/C/D"); + } + + #[test] + fn block_decode_rate_osd1_vs_osd2() { + // Measure how many blocks with exactly 2 bit errors are recovered + // by OSD(2) vs the number that would succeed at OSD(1) (none for 2-bit errors). + // + // For a valid block with 2 bit errors at the two *least* confident + // positions, OSD(2) should always recover it; OSD(1) should never. + let offsets = [OFFSET_A, OFFSET_B, OFFSET_C, OFFSET_CP, OFFSET_D]; + let data_values: [u16; 5] = [0x1111, 0x2222, 0x4865, 0xBEEF, 0xCAFE]; + + let mut osd1_ok = 0u32; + let mut osd2_ok = 0u32; + let total = offsets.len() * data_values.len(); + + for &offset in &offsets { + for &data in &data_values { + let word = encode_block(data, offset); + // Flip bits k=0 and k=25 (spread across the word). + let corrupted = word ^ (1 << 25) ^ (1 << 0); + let mut soft = [1.0f32; 26]; + soft[0] = 0.01; // very uncertain + soft[25] = 0.01; + + // OSD(1): should fail for 2-bit errors. + let osd1_result = { + if decode_block(corrupted).is_some() { + Some(()) // d0 hit (unexpected but count it) + } else { + (0..26usize).find_map(|k| decode_block(corrupted ^ (1 << (25 - k)))).map(|_| ()) + } + }; + if osd1_result.is_some() { + osd1_ok += 1; + } + + // OSD(2). + if decode_block_soft(corrupted, &soft).is_some() { + osd2_ok += 1; + } + } + } + + // OSD(2) must recover at least 80% of cleanly 2-bit-corrupted blocks. + assert!( + osd2_ok >= (total as u32 * 8 / 10), + "OSD(2) recovery rate = {}/{total} (< 80%)", + osd2_ok + ); + // OSD(1) should not recover any (all are genuine 2-bit errors). + assert_eq!(osd1_ok, 0, "OSD(1) should not recover 2-bit errors"); + } + + // ----------------------------------------------------------------------- + // Costas loop convergence + // ----------------------------------------------------------------------- + + #[test] + fn costas_tracks_without_diverging_on_clean_signal() { + // Feed a clean RDS signal through the decoder (no pilot ref) and + // verify that at least some PI data is recovered, proving the Costas + // loop stays coherent rather than losing lock permanently. + let sample_rate = 240_000.0f32; + let pi = 0x7777u16; + + let mut words: Vec = Vec::new(); + for seg in 0..4u8 { + let g = group_0a(pi, seg, [b'C' + seg, b'D' + seg], 5); + words.extend_from_slice(&g); + } + // 60× repetitions to give Costas plenty of time to acquire. + let words: Vec = words.iter().copied().cycle().take(words.len() * 60).collect(); + + let chips = blocks_to_chips(&words); + let signal = chips_to_rds_signal(&chips, sample_rate); + + let mut dec = RdsDecoder::new(sample_rate as u32); + let mut pi_correct = 0u32; + let mut pi_total = 0u32; + for &s in &signal { + if let Some(state) = dec.process_sample(s, 1.0) { + if state.pi.is_some() { + pi_total += 1; + if state.pi == Some(pi) { + pi_correct += 1; + } + } + } + } + assert!( + pi_correct > 0, + "Costas should converge and produce correct PI (got {pi_correct}/{pi_total})" + ); + } + + // ----------------------------------------------------------------------- + // Noise rejection + // ----------------------------------------------------------------------- + + #[test] + fn pure_noise_produces_zero_pi_decodes() { + // Feed 0.5 seconds of white noise (no RDS signal) through the decoder. + // The decoder must not report any PI (false positive). + // + // Note: with OSD(2) active in locked mode, the lock gate requires + // Block A to be acquired first (hard or OSD-1 decode in search mode), + // which keeps the false-acquisition rate low even at OSD(2). + let sample_rate = 240_000.0f32; + let n_samples = (sample_rate * 0.5) as usize; + let mut rng = 0xFEED_FACE_DEAD_BEEFu64; + let mut noise: Vec = (0..n_samples).map(|_| gaussian(&mut rng)).collect(); + // Scale noise to unit power. + let pwr = noise.iter().map(|x| x * x).sum::() / n_samples as f32; + let scale = pwr.sqrt().recip(); + noise.iter_mut().for_each(|x| *x *= scale); + + let mut dec = RdsDecoder::new(sample_rate as u32); + let mut false_pi = 0u32; + for &s in &noise { + if let Some(state) = dec.process_sample(s, 1.0) { + if state.pi.is_some() { + false_pi += 1; + } + } + } + assert_eq!( + false_pi, 0, + "Pure noise generated {false_pi} false PI reports" + ); + } + + // ----------------------------------------------------------------------- + // PI accumulation + // ----------------------------------------------------------------------- + + #[test] + fn pi_accumulation_corrects_weak_pi_after_threshold() { + // The PI LLR accumulator (Tech 6) should vote out a one-bit error + // in the PI field after PI_ACC_THRESHOLD observations. + let real_pi: u16 = 0x9420; + let bad_pi: u16 = real_pi ^ 0x0001; // one bit wrong in LSB + + let pty = 10u8; + let mut cand = Candidate::new(240_000.0, 0.0); + + // Send PI_ACC_THRESHOLD groups; each Block A carries the correct PI + // but with a very low soft confidence on the corrupted bit position. + for i in 0..(PI_ACC_THRESHOLD + 1) { + let block_a = encode_block(real_pi, OFFSET_A); + let block_b = encode_block(u16::from(pty) << 5, OFFSET_B); + let block_c = encode_block(0, OFFSET_C); + let block_d = encode_block(u16::from_be_bytes(*b"OK"), OFFSET_D); + + for bit_idx in (0..26).rev() { + let bit = ((block_a >> bit_idx) & 1) as u8; + // Bit 0 (LSB of PI) has low confidence. + let conf = if bit_idx == 0 { 0.1 } else { 1.0 }; + let _ = cand.push_bit_soft(bit, conf); + } + for bit_idx in (0..26).rev() { + let _ = cand.push_bit_soft(((block_b >> bit_idx) & 1) as u8, 1.0); + } + for bit_idx in (0..26).rev() { + let _ = cand.push_bit_soft(((block_c >> bit_idx) & 1) as u8, 1.0); + } + let mut last = None; + for bit_idx in (0..26).rev() { + last = cand.push_bit_soft(((block_d >> bit_idx) & 1) as u8, 1.0); + } + let _ = (i, last, bad_pi); // silence unused warnings + } + + // After threshold groups, the accumulated PI should converge to real_pi. + let pi = cand.state.pi.expect("PI should be set after accumulation"); + assert_eq!( + pi, real_pi, + "Accumulated PI {pi:#06x} should converge to {real_pi:#06x}" + ); + } } diff --git a/src/trx-server/trx-backend/trx-backend-soapysdr/src/demod/wfm.rs b/src/trx-server/trx-backend/trx-backend-soapysdr/src/demod/wfm.rs index 60d637c..fd822d9 100644 --- a/src/trx-server/trx-backend/trx-backend-soapysdr/src/demod/wfm.rs +++ b/src/trx-server/trx-backend/trx-backend-soapysdr/src/demod/wfm.rs @@ -10,7 +10,13 @@ use super::{math::demod_fm_with_prev, DcBlocker}; const RDS_SUBCARRIER_HZ: f32 = 57_000.0; /// Tech 2: pilot lock level above which the ×3 pilot reference is used. -const PILOT_LOCK_THRESHOLD: f32 = 0.25; +/// Effective pilot coherence threshold ≈ ONSET + THRESHOLD × 0.2 = 0.36. +const PILOT_LOCK_THRESHOLD: f32 = 0.20; +/// Coherence below which pilot_lock contribution is zero (linear ramp 0→1 +/// over the range [ONSET, ONSET+0.2]). Lower value → pilot ref used on +/// weaker stations; risk: noisier reference. 0.30 vs original 0.40 means +/// we engage at coherence ≥ 0.36 instead of ≥ 0.45. +const PILOT_LOCK_ONSET: f32 = 0.30; /// Tech 9: number of complex CMA equalizer taps. const CMA_N_TAPS: usize = 8; /// Tech 9: CMA LMS step size. @@ -722,7 +728,7 @@ impl WfmStereoDecoder { let avg_mag = self.detect_pilot_mag_acc * inv_n; let avg_abs = self.detect_pilot_abs_acc * inv_n; let pilot_coherence = (avg_mag / (avg_abs + 1e-4)).clamp(0.0, 1.0); - let pilot_lock = ((pilot_coherence - 0.4) / 0.2).clamp(0.0, 1.0); + let pilot_lock = ((pilot_coherence - PILOT_LOCK_ONSET) / 0.2).clamp(0.0, 1.0); self.pilot_lock_level += 0.12 * (pilot_lock - self.pilot_lock_level); let stereo_drive = (avg_mag * pilot_lock * 120.0).clamp(0.0, 1.0); let detect_coeff = if stereo_drive > self.stereo_detect_level {