From 3877de3c4f8d63991e304176863e31fedf0c7a08 Mon Sep 17 00:00:00 2001 From: Stan Grams Date: Fri, 20 Mar 2026 08:16:33 +0100 Subject: [PATCH] [fix](trx-wspr): rewrite decoder with soft-decision Fano from WSJT-X reference Replace broken hard-decision pipeline with proper soft-decision decoding matching the WSJT-X wsprd reference implementation. Key fixes: soft-symbol demodulation using amplitude differences, soft-decision Fano decoder with Es/No=6dB metric table and delta=60 threshold, deinterleave preserving soft values instead of extracting hard bits, convolutional tail constraint, and normalized sync correlation scoring. Co-Authored-By: Claude Opus 4.6 Signed-off-by: Stan Grams --- src/decoders/trx-wspr/src/decoder.rs | 226 ++++++++++----- src/decoders/trx-wspr/src/protocol.rs | 397 +++++++++++++++++++++----- 2 files changed, 480 insertions(+), 143 deletions(-) diff --git a/src/decoders/trx-wspr/src/decoder.rs b/src/decoders/trx-wspr/src/decoder.rs index 6bfca16..ea0df63 100644 --- a/src/decoders/trx-wspr/src/decoder.rs +++ b/src/decoders/trx-wspr/src/decoder.rs @@ -18,22 +18,24 @@ const BASE_SEARCH_MAX_HZ: f32 = 1800.0; const BASE_SEARCH_STEP_HZ: f32 = 2.0; const FINE_SEARCH_STEP_HZ: f32 = 0.25; -// Timing offset search: search ±2s in 0.5s steps (4800 samples at 12 kHz) +// Timing offset search: search ±2s in 0.5s steps (6000 samples at 12 kHz) const DT_SEARCH_RANGE_SAMPLES: isize = 2 * WSPR_SAMPLE_RATE as isize; const DT_SEARCH_STEP_SAMPLES: isize = (WSPR_SAMPLE_RATE as isize) / 2; // Number of top frequency candidates to try full decode on const MAX_FREQ_CANDIDATES: usize = 8; -// Minimum sync correlation score to attempt a full decode. Candidates below -// this threshold are almost certainly noise and skipping them avoids expensive -// Fano decode attempts that would produce false positives. -const MIN_SYNC_SCORE: f32 = 10.0; +// Minimum normalized sync correlation score to attempt decode. +// Matches reference wsprd minsync1=0.10. +const MIN_SYNC_SCORE: f32 = 0.10; + +// Soft-symbol normalization factor (reference wsprd: symfac=50) +const SYMFAC: f32 = 50.0; /// WSPR sync vector (162 bits). symbol = sync[i] + 2*data[i]. /// The LSB of each received symbol should match this pattern. #[rustfmt::skip] -const SYNC_VECTOR: [u8; 162] = [ +pub(crate) const SYNC_VECTOR: [u8; 162] = [ 1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,0,0,0, 0,0,0,0,1,0,0,1,0,1,0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1, 1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,1,1,0,0,0,1, @@ -54,6 +56,11 @@ pub struct WsprDecoder { min_rms: f32, } +struct DemodOutput { + soft_symbols: [u8; WSPR_SYMBOL_COUNT], + snr_db: f32, +} + impl WsprDecoder { pub fn new() -> Result { Ok(Self { min_rms: 0.0005 }) @@ -128,15 +135,18 @@ impl WsprDecoder { let mut results = Vec::new(); let mut seen_messages = std::collections::HashSet::new(); - for &(freq, dt_samples, score) in candidates.iter().take(MAX_FREQ_CANDIDATES) { - if score < MIN_SYNC_SCORE { - break; // candidates are sorted by score, no point continuing - } + for &(freq, dt_samples, _score) in candidates.iter().take(MAX_FREQ_CANDIDATES) { let start = (EXPECTED_SIGNAL_START_SAMPLES as isize + dt_samples) as usize; let signal = &samples[start..start + WSPR_SIGNAL_SAMPLES]; - let demod = demodulate_symbols(signal, freq); - if let Some(decoded) = protocol::decode_symbols(&demod.symbols) { + // Use normalized sync score for threshold check + let norm_score = sync_correlation_score_normalized(signal, freq); + if norm_score < MIN_SYNC_SCORE { + continue; + } + + let demod = demodulate_soft_symbols(signal, freq); + if let Some(decoded) = protocol::decode_symbols(&demod.soft_symbols) { if seen_messages.insert(decoded.message.clone()) { let dt_s = dt_samples as f32 / WSPR_SAMPLE_RATE as f32; results.push(WsprDecodeResult { @@ -153,49 +163,81 @@ impl WsprDecoder { } } -#[derive(Debug, Clone)] -struct DemodOutput { - symbols: Vec, - snr_db: f32, -} - -/// Score a candidate base frequency by correlating detected symbol LSBs with -/// the known WSPR sync vector. Higher score = better match. +/// Score a candidate base frequency by correlating detected tone amplitudes +/// with the known WSPR sync vector. Uses amplitude (sqrt of power) and +/// normalizes by total power, matching the reference wsprd implementation. +/// Higher score = better match. Range approximately [0.0, 1.0]. fn sync_correlation_score(signal: &[f32], base_hz: f32) -> f32 { let nsyms = WSPR_SYMBOL_COUNT.min(signal.len() / WSPR_SYMBOL_SAMPLES); - let mut score = 0.0_f32; + let mut ss = 0.0_f32; + let sr = WSPR_SAMPLE_RATE as f32; + for (sym, &sync_bit) in SYNC_VECTOR.iter().enumerate().take(nsyms) { let off = sym * WSPR_SYMBOL_SAMPLES; let frame = &signal[off..off + WSPR_SYMBOL_SAMPLES]; - // Sum power in even tones (0,2) vs odd tones (1,3) - let p0 = goertzel_power(frame, base_hz, WSPR_SAMPLE_RATE as f32); - let p2 = goertzel_power( - frame, - base_hz + 2.0 * TONE_SPACING_HZ, - WSPR_SAMPLE_RATE as f32, - ); - let p1 = goertzel_power(frame, base_hz + TONE_SPACING_HZ, WSPR_SAMPLE_RATE as f32); - let p3 = goertzel_power( - frame, - base_hz + 3.0 * TONE_SPACING_HZ, - WSPR_SAMPLE_RATE as f32, - ); - let even_power = p0 + p2; // tones with LSB=0 - let odd_power = p1 + p3; // tones with LSB=1 + // Compute amplitude (sqrt of power) at each of the 4 FSK tones + let p0 = goertzel_power(frame, base_hz, sr).sqrt(); + let p1 = goertzel_power(frame, base_hz + TONE_SPACING_HZ, sr).sqrt(); + let p2 = goertzel_power(frame, base_hz + 2.0 * TONE_SPACING_HZ, sr).sqrt(); + let p3 = goertzel_power(frame, base_hz + 3.0 * TONE_SPACING_HZ, sr).sqrt(); - // Correlate with sync vector: sync=1 means odd tone expected + // Correlate with sync vector: (p1+p3)-(p0+p2) weighted by (2*sync-1) + let cmet = (p1 + p3) - (p0 + p2); if sync_bit == 1 { - score += odd_power - even_power; + ss += cmet; } else { - score += even_power - odd_power; + ss -= cmet; } } - score + + // Raw (unnormalized) score for candidate ranking. At frequencies with no + // signal, amplitude differences are near zero so raw score is naturally low. + // Normalized threshold check is applied separately in decode_slot. + ss } -fn demodulate_symbols(signal: &[f32], base_hz: f32) -> DemodOutput { - let mut symbols = Vec::with_capacity(WSPR_SYMBOL_COUNT); +/// Compute the normalized sync score (ss/totp) for threshold comparison. +fn sync_correlation_score_normalized(signal: &[f32], base_hz: f32) -> f32 { + let nsyms = WSPR_SYMBOL_COUNT.min(signal.len() / WSPR_SYMBOL_SAMPLES); + let mut ss = 0.0_f32; + let mut totp = 0.0_f32; + let sr = WSPR_SAMPLE_RATE as f32; + + for (sym, &sync_bit) in SYNC_VECTOR.iter().enumerate().take(nsyms) { + let off = sym * WSPR_SYMBOL_SAMPLES; + let frame = &signal[off..off + WSPR_SYMBOL_SAMPLES]; + + let p0 = goertzel_power(frame, base_hz, sr).sqrt(); + let p1 = goertzel_power(frame, base_hz + TONE_SPACING_HZ, sr).sqrt(); + let p2 = goertzel_power(frame, base_hz + 2.0 * TONE_SPACING_HZ, sr).sqrt(); + let p3 = goertzel_power(frame, base_hz + 3.0 * TONE_SPACING_HZ, sr).sqrt(); + + let cmet = (p1 + p3) - (p0 + p2); + if sync_bit == 1 { + ss += cmet; + } else { + ss -= cmet; + } + totp += p0 + p1 + p2 + p3; + } + + if totp > 0.0 { + ss / totp + } else { + 0.0 + } +} + +/// Produce soft-decision symbols from a signal slice. +/// +/// Each soft symbol is an unsigned byte (0-255) where 128 = no confidence, +/// values above 128 mean data bit is likely 1, below 128 means likely 0. +/// +/// This matches the reference wsprd `sync_and_demodulate` mode=2 output. +fn demodulate_soft_symbols(signal: &[f32], base_hz: f32) -> DemodOutput { + let sr = WSPR_SAMPLE_RATE as f32; + let mut fsymb = [0.0_f32; WSPR_SYMBOL_COUNT]; let mut signal_sum = 0.0_f32; let mut noise_sum = 0.0_f32; @@ -203,44 +245,56 @@ fn demodulate_symbols(signal: &[f32], base_hz: f32) -> DemodOutput { let off = sym * WSPR_SYMBOL_SAMPLES; let frame = &signal[off..off + WSPR_SYMBOL_SAMPLES]; - let mut tone_power = [0.0_f32; 4]; - for (i, power) in tone_power.iter_mut().enumerate() { - let hz = base_hz + i as f32 * TONE_SPACING_HZ; - *power = goertzel_power(frame, hz, WSPR_SAMPLE_RATE as f32); + // Compute amplitude (sqrt of power) at each tone — matches reference + let p0 = goertzel_power(frame, base_hz, sr).sqrt(); + let p1 = goertzel_power(frame, base_hz + TONE_SPACING_HZ, sr).sqrt(); + let p2 = goertzel_power(frame, base_hz + 2.0 * TONE_SPACING_HZ, sr).sqrt(); + let p3 = goertzel_power(frame, base_hz + 3.0 * TONE_SPACING_HZ, sr).sqrt(); + + // Soft metric for the data bit: + // sync=1 → data bit selects tone 1 (data=0) vs tone 3 (data=1) + // sync=0 → data bit selects tone 0 (data=0) vs tone 2 (data=1) + // Positive fsymb means data_bit=1 is more likely. + if SYNC_VECTOR[sym] == 1 { + fsymb[sym] = p3 - p1; + } else { + fsymb[sym] = p2 - p0; } - let mut best_idx = 0_u8; - let mut best_pow = tone_power[0]; - for (idx, p) in tone_power.iter().enumerate().skip(1) { - if *p > best_pow { - best_pow = *p; - best_idx = idx as u8; - } - } + // SNR estimation: signal = best tone power, noise = out-of-band + let best_amp = p0.max(p1).max(p2).max(p3); + signal_sum += best_amp * best_amp; - symbols.push(best_idx); - signal_sum += best_pow; - - let noise_a = goertzel_power( - frame, - base_hz - 8.0 * TONE_SPACING_HZ, - WSPR_SAMPLE_RATE as f32, - ); - let noise_b = goertzel_power( - frame, - base_hz + 12.0 * TONE_SPACING_HZ, - WSPR_SAMPLE_RATE as f32, - ); + let noise_a = goertzel_power(frame, base_hz - 8.0 * TONE_SPACING_HZ, sr); + let noise_b = goertzel_power(frame, base_hz + 12.0 * TONE_SPACING_HZ, sr); noise_sum += (noise_a + noise_b) * 0.5; } - let signal_avg = signal_sum / WSPR_SYMBOL_COUNT as f32; - let noise_avg = (noise_sum / WSPR_SYMBOL_COUNT as f32).max(1e-12); + // Normalize: zero-mean, scale by symfac/stddev, clip to [-128,127], bias to [0,255] + let n = WSPR_SYMBOL_COUNT as f32; + let mean = fsymb.iter().sum::() / n; + let var = fsymb.iter().map(|&x| (x - mean) * (x - mean)).sum::() / n; + let fac = var.sqrt().max(1e-12); + + let mut soft_symbols = [128u8; WSPR_SYMBOL_COUNT]; + for i in 0..WSPR_SYMBOL_COUNT { + let v = SYMFAC * fsymb[i] / fac; + let v = v.clamp(-128.0, 127.0); + soft_symbols[i] = (v + 128.0) as u8; + } + + // SNR estimate + let signal_avg = signal_sum / n; + let noise_avg = (noise_sum / n).max(1e-12); let snr_db = 10.0 * (signal_avg / noise_avg).max(1e-12).log10(); - DemodOutput { symbols, snr_db } + DemodOutput { + soft_symbols, + snr_db, + } } +/// Goertzel algorithm: compute power at a specific frequency in a windowed frame. fn goertzel_power(frame: &[f32], target_hz: f32, sample_rate: f32) -> f32 { let n = frame.len() as f32; let k = (0.5 + (n * target_hz / sample_rate)).floor(); @@ -373,4 +427,36 @@ mod tests { "correct={correct_score}, wrong={wrong_score}" ); } + + #[test] + fn normalized_sync_score_is_bounded() { + let base_hz = 1500.0_f32; + + // Generate a perfect synthetic WSPR signal + let mut signal = vec![0.0_f32; WSPR_SIGNAL_SAMPLES]; + for (sym, sync_tone) in SYNC_VECTOR + .iter() + .copied() + .enumerate() + .take(WSPR_SYMBOL_COUNT) + { + // Use sync_tone as the only varying bit to maximize sync metric + let freq = base_hz + sync_tone as f32 * TONE_SPACING_HZ; + let begin = sym * WSPR_SYMBOL_SAMPLES; + for i in 0..WSPR_SYMBOL_SAMPLES { + let t = i as f32 / WSPR_SAMPLE_RATE as f32; + signal[begin + i] = (2.0 * std::f32::consts::PI * freq * t).sin() * 0.2; + } + } + + let score = sync_correlation_score_normalized(&signal, base_hz); + // Normalized score should be positive and bounded + assert!(score > 0.0, "score should be positive: {score}"); + assert!(score <= 1.0, "score should be <= 1.0: {score}"); + // For a signal where only sync tones are present, score should be high + assert!( + score > MIN_SYNC_SCORE, + "score {score} should exceed threshold {MIN_SYNC_SCORE}" + ); + } } diff --git a/src/decoders/trx-wspr/src/protocol.rs b/src/decoders/trx-wspr/src/protocol.rs index 2095468..ce6d26d 100644 --- a/src/decoders/trx-wspr/src/protocol.rs +++ b/src/decoders/trx-wspr/src/protocol.rs @@ -13,6 +13,63 @@ const POLY2: u32 = 0xE4613C47; const NBITS: usize = 81; // 50 payload bits + 31 convolutional flush bits const NSYMS: usize = 162; +// Fano decoder parameters (matching reference wsprd) +const FANO_DELTA: i32 = 60; +const FANO_MAX_CYCLES_PER_BIT: usize = 10_000; +const FANO_BIAS: f32 = 0.45; + +/// Soft-decision metric table for the Fano decoder. +/// Es/No = 6 dB log-likelihood ratio table from WSJT-X reference (metric_tables[2]). +#[allow(clippy::approx_constant)] +#[rustfmt::skip] +const METRIC_TABLE: [f32; 256] = [ + 0.9999, 0.9998, 0.9998, 0.9998, 0.9998, 0.9998, 0.9997, 0.9997, + 0.9997, 0.9997, 0.9997, 0.9996, 0.9996, 0.9996, 0.9995, 0.9995, + 0.9994, 0.9994, 0.9994, 0.9993, 0.9993, 0.9992, 0.9991, 0.9991, + 0.9990, 0.9989, 0.9988, 0.9988, 0.9988, 0.9986, 0.9985, 0.9984, + 0.9983, 0.9982, 0.9980, 0.9979, 0.9977, 0.9976, 0.9974, 0.9971, + 0.9969, 0.9968, 0.9965, 0.9962, 0.9960, 0.9957, 0.9953, 0.9950, + 0.9947, 0.9941, 0.9937, 0.9933, 0.9928, 0.9922, 0.9917, 0.9911, + 0.9904, 0.9897, 0.9890, 0.9882, 0.9874, 0.9863, 0.9855, 0.9843, + 0.9832, 0.9819, 0.9806, 0.9792, 0.9777, 0.9760, 0.9743, 0.9724, + 0.9704, 0.9683, 0.9659, 0.9634, 0.9609, 0.9581, 0.9550, 0.9516, + 0.9481, 0.9446, 0.9406, 0.9363, 0.9317, 0.9270, 0.9218, 0.9160, + 0.9103, 0.9038, 0.8972, 0.8898, 0.8822, 0.8739, 0.8647, 0.8554, + 0.8457, 0.8357, 0.8231, 0.8115, 0.7984, 0.7854, 0.7704, 0.7556, + 0.7391, 0.7210, 0.7038, 0.6840, 0.6633, 0.6408, 0.6174, 0.5939, + 0.5678, 0.5410, 0.5137, 0.4836, 0.4524, 0.4193, 0.3850, 0.3482, + 0.3132, 0.2733, 0.2315, 0.1891, 0.1435, 0.0980, 0.0493, 0.0000, + -0.0510, -0.1052, -0.1593, -0.2177, -0.2759, -0.3374, -0.4005, -0.4599, + -0.5266, -0.5935, -0.6626, -0.7328, -0.8051, -0.8757, -0.9498, -1.0271, + -1.1019, -1.1816, -1.2642, -1.3459, -1.4295, -1.5077, -1.5958, -1.6818, + -1.7647, -1.8548, -1.9387, -2.0295, -2.1152, -2.2154, -2.3011, -2.3904, + -2.4820, -2.5786, -2.6730, -2.7652, -2.8616, -2.9546, -3.0526, -3.1445, + -3.2445, -3.3416, -3.4357, -3.5325, -3.6324, -3.7313, -3.8225, -3.9209, + -4.0248, -4.1278, -4.2261, -4.3193, -4.4220, -4.5262, -4.6214, -4.7242, + -4.8234, -4.9245, -5.0298, -5.1250, -5.2232, -5.3267, -5.4332, -5.5342, + -5.6431, -5.7270, -5.8401, -5.9350, -6.0407, -6.1418, -6.2363, -6.3384, + -6.4536, -6.5429, -6.6582, -6.7433, -6.8438, -6.9478, -7.0789, -7.1894, + -7.2714, -7.3815, -7.4810, -7.5575, -7.6852, -7.8071, -7.8580, -7.9724, + -8.1000, -8.2207, -8.2867, -8.4017, -8.5287, -8.6347, -8.7082, -8.8319, + -8.9448, -9.0355, -9.1885, -9.2095, -9.2863, -9.4186, -9.5064, -9.6386, + -9.7207, -9.8286, -9.9453,-10.0701,-10.1735,-10.3001,-10.2858,-10.5427, + -10.5982,-10.7361,-10.7042,-10.9212,-11.0097,-11.0469,-11.1155,-11.2812, + -11.3472,-11.4988,-11.5327,-11.6692,-11.9376,-11.8606,-12.1372,-13.2539, +]; + +/// Build the integer metric table for the soft-decision Fano decoder. +/// +/// `mettab[0][rx]` = metric when expected coded bit is 0, received symbol is `rx` +/// `mettab[1][rx]` = metric when expected coded bit is 1, received symbol is `rx` +fn build_mettab() -> [[i32; 256]; 2] { + let mut mettab = [[0i32; 256]; 2]; + for i in 0..256 { + mettab[0][i] = (10.0 * (METRIC_TABLE[i] - FANO_BIAS)).round() as i32; + mettab[1][i] = (10.0 * (METRIC_TABLE[255 - i] - FANO_BIAS)).round() as i32; + } + mettab +} + /// Reverse the bits of an 8-bit value. fn rev8(mut b: u8) -> u8 { let mut r = 0u8; @@ -23,99 +80,159 @@ fn rev8(mut b: u8) -> u8 { r } -/// Extract hard data bits from 4-FSK symbols then deinterleave. +/// Deinterleave soft symbols by permuting their order via bit-reversal of indices. /// -/// In WSPR: symbol = sync_bit + 2*data_bit, so data_bit = symbol >> 1. -/// The 162 data bits are reordered via bit-reversal of 8-bit indices. -/// -/// The interleaving places coded bit p at transmitted position j = rev8(i) -/// (for each i in 0..255 where rev8(i) < 162). Deinterleaving reverses -/// this: coded[p] = transmitted[j] >> 1. +/// Unlike the old hard-decision version, this does NOT extract data bits — the +/// soft values (0-255, centered at 128) are preserved as-is. The Fano decoder +/// interprets them directly via the metric table. fn deinterleave(symbols: &[u8]) -> [u8; NSYMS] { - let mut out = [0u8; NSYMS]; + let mut out = [128u8; NSYMS]; // default to "no confidence" let mut p = 0usize; - for i in 0u8..=255 { - let j = rev8(i) as usize; + for i in 0u16..=255 { + let j = rev8(i as u8) as usize; if j < NSYMS { - out[p] = symbols[j] >> 1; + out[p] = if j < symbols.len() { symbols[j] } else { 128 }; p += 1; } } out } -/// Fano sequential decoder for the K=32, rate-1/2 convolutional code. +/// Compute the 2-bit convolutional encoder output for a given encoder state. /// -/// `coded[2*k]` and `coded[2*k+1]` are the two received bits for input bit k. -/// Returns the 81 decoded bits (first 50 are payload), or None if it cannot -/// converge within the iteration budget. -fn fano_decode(coded: &[u8; NSYMS]) -> Option<[u8; NBITS]> { - const MAX_CYCLES: usize = 100_000; +/// Returns a value 0-3 where: +/// bit 1 (2's place) = parity(state & POLY1) +/// bit 0 (1's place) = parity(state & POLY2) +fn encode_sym(state: u32) -> u32 { + let p1 = (state & POLY1).count_ones() & 1; + let p2 = (state & POLY2).count_ones() & 1; + (p1 << 1) | p2 +} - // Hard-decision branch metric: +1 per matching bit, -1 per mismatch. - let bm = |pos: usize, state: u32, bit: u32| -> i32 { - let ns = (state << 1) | bit; - let c0 = ((ns & POLY1).count_ones() & 1) as u8; - let c1 = ((ns & POLY2).count_ones() & 1) as u8; - let r0 = coded[2 * pos]; - let r1 = coded[2 * pos + 1]; - (if c0 == r0 { 1 } else { -1 }) + (if c1 == r1 { 1 } else { -1 }) - }; +/// Soft-decision Fano sequential decoder for K=32, rate-1/2 convolutional code. +/// +/// Closely follows the reference implementation from WSJT-X (fano.c by Phil Karn, KA9Q). +/// +/// Input: 162 deinterleaved soft-decision symbols (0-255, 128=no confidence). +/// Symbols are read in pairs: `symbols[2k]` and `symbols[2k+1]` are the two +/// coded bits for input bit k. +/// +/// Output: 81 decoded bits (first 50 are payload), or None on timeout. +fn fano_decode(symbols: &[u8; NSYMS]) -> Option<[u8; NBITS]> { + let mettab = build_mettab(); + let max_cycles = FANO_MAX_CYCLES_PER_BIT * NBITS; + let tail_start = NBITS - 31; // position 50: first tail bit - let mut bits = [0u8; NBITS]; - let mut states = [0u32; NBITS + 1]; - let mut metrics = [0i32; NBITS + 1]; - // 0 = not yet visited, 1 = tried best branch, 2 = tried both branches - let mut visit = [0u8; NBITS]; + // Precompute all 4 branch metrics for each bit position. + // metrics[k][sym_pair] where sym_pair encodes (expected_bit0, expected_bit1): + // 0 = (0,0), 1 = (0,1), 2 = (1,0), 3 = (1,1) + let mut metrics = [[0i32; 4]; NBITS]; + for k in 0..NBITS { + let s0 = symbols[2 * k] as usize; + let s1 = symbols[2 * k + 1] as usize; + metrics[k][0] = mettab[0][s0] + mettab[0][s1]; + metrics[k][1] = mettab[0][s0] + mettab[1][s1]; + metrics[k][2] = mettab[1][s0] + mettab[0][s1]; + metrics[k][3] = mettab[1][s0] + mettab[1][s1]; + } + + // Per-node state + let mut encstate = [0u32; NBITS + 1]; + let mut gamma = [0i64; NBITS + 1]; // cumulative path metric + let mut tm = [[0i32; 2]; NBITS]; // sorted branch metrics [best, second] + let mut branch_i = [0u8; NBITS]; // 0 = trying best branch, 1 = trying second let mut pos: usize = 0; - let mut threshold: i32 = 0; + let mut t: i64 = 0; // threshold - for _ in 0..MAX_CYCLES { - if pos == NBITS { + // Initialize root node: compute and sort branch metrics + let lsym = encode_sym(encstate[0]) as usize; + let m0 = metrics[0][lsym]; + let m1 = metrics[0][3 ^ lsym]; + if m0 > m1 { + tm[0] = [m0, m1]; + } else { + tm[0] = [m1, m0]; + encstate[0] |= 1; // 1-branch is better; encode choice in LSB + } + branch_i[0] = 0; + + for _cycle in 0..max_cycles { + if pos >= NBITS { break; } - let state = states[pos]; - let base = metrics[pos]; - let m0 = base + bm(pos, state, 0); - let m1 = base + bm(pos, state, 1); - - // b0/mv0 = best branch; b1/mv1 = second-best - let (b0, mv0, b1, mv1) = if m0 >= m1 { - (0u8, m0, 1u8, m1) - } else { - (1u8, m1, 0u8, m0) - }; - - // Which branch to try at this visit count? - let chosen = match visit[pos] { - 0 if mv0 >= threshold => Some((b0, mv0)), - 1 if mv1 >= threshold => Some((b1, mv1)), - _ => None, - }; - - if let Some((b, m)) = chosen { - visit[pos] = if visit[pos] == 0 { 1 } else { 2 }; - bits[pos] = b; - states[pos + 1] = (state << 1) | b as u32; - metrics[pos + 1] = m; - pos += 1; - } else { - // Both branches at this position fail the threshold: backtrack. - visit[pos] = 0; - if pos == 0 { - threshold -= 1; - } else { - pos -= 1; - // Parent node (visit[pos] == 1 or 2) will try the next branch - // on the following iteration, or backtrack further. + // Look forward: try current branch + let ngamma = gamma[pos] + tm[pos][branch_i[pos] as usize] as i64; + if ngamma >= t { + // Acceptable — tighten threshold if this is a first visit + if gamma[pos] < t + FANO_DELTA as i64 { + while ngamma >= t + FANO_DELTA as i64 { + t += FANO_DELTA as i64; + } } + + // Move forward + gamma[pos + 1] = ngamma; + encstate[pos + 1] = encstate[pos] << 1; + pos += 1; + + if pos >= NBITS { + break; // Done! + } + + // Compute and sort metrics at the new position + let lsym = encode_sym(encstate[pos]) as usize; + if pos >= tail_start { + // Tail must be all zeros — only consider 0-branch + tm[pos] = [metrics[pos][lsym], i32::MIN]; + } else { + let m0 = metrics[pos][lsym]; + let m1 = metrics[pos][3 ^ lsym]; + if m0 > m1 { + tm[pos] = [m0, m1]; + } else { + tm[pos] = [m1, m0]; + encstate[pos] |= 1; // mark 1-branch as better + } + } + branch_i[pos] = 0; + continue; + } + + // Threshold violated — look backward + loop { + if pos == 0 || gamma[pos - 1] < t { + // Can't back up (at root, or parent's metric below threshold). + // Relax threshold and reset to best branch at current position. + t -= FANO_DELTA as i64; + if branch_i[pos] != 0 { + branch_i[pos] = 0; + encstate[pos] ^= 1; + } + break; + } + // Back up to parent + pos -= 1; + if pos < tail_start && branch_i[pos] != 1 { + // Try second branch at this position + branch_i[pos] = 1; + encstate[pos] ^= 1; + break; + } + // Already tried both branches (or in tail) — keep backing up } } if pos < NBITS { - return None; + return None; // Timeout + } + + // Extract decoded bits from encoder states. + // At each position k, the LSB of encstate[k] is the chosen input bit. + let mut bits = [0u8; NBITS]; + for k in 0..NBITS { + bits[k] = (encstate[k] & 1) as u8; } Some(bits) } @@ -226,7 +343,12 @@ fn unpack_message(bits: &[u8; NBITS]) -> Option { Some(format!("{} {} {}", callsign, grid, power_dbm)) } -/// Attempt protocol-level decode from 162 4-FSK symbols. +/// Attempt protocol-level decode from 162 soft-decision symbols. +/// +/// Input: 162 bytes where each value is a soft-decision symbol (0-255): +/// 0 = high confidence that data bit is 0 +/// 128 = no confidence +/// 255 = high confidence that data bit is 1 pub fn decode_symbols(symbols: &[u8]) -> Option { if symbols.len() < NSYMS { return None; @@ -240,6 +362,7 @@ pub fn decode_symbols(symbols: &[u8]) -> Option { #[cfg(test)] mod tests { use super::*; + use crate::decoder::SYNC_VECTOR; /// Encode a WSPR callsign+grid+power into N1/M1/power_code, then round-trip /// through `unpack_message` to verify the pack/unpack formulas are inverse. @@ -254,9 +377,6 @@ mod tests { let idx37 = |c: u8| cs37.iter().position(|&x| x == c).unwrap() as u32; let idx27 = |c: u8| cs27.iter().position(|&x| x == c).unwrap() as u32; - // "K1JT " → c0='K'=21, c1='1'=2, c2='J'-no, wait: position 2 must be digit - // Standard WSPR normalises so that position 2 is the digit. - // "K1JT " has digit at position 1, not 2 → needs prefix space: " K1JT " // " K1JT ": c0=' '=0, c1='K'=21, c2='1', c3='J'=10, c4='T'=20, c5=' '=0 let c0 = idx37(b' '); let c1 = idx37(b'K'); @@ -295,4 +415,135 @@ mod tests { assert!(msg.contains("FN20"), "grid not found in '{}'", msg); assert!(msg.contains("37"), "power not found in '{}'", msg); } + + /// Convolutionally encode 81 bits → 162 coded bits (for testing). + fn convolutional_encode(input: &[u8; NBITS]) -> [u8; NSYMS] { + let mut coded = [0u8; NSYMS]; + let mut encstate: u32 = 0; + for k in 0..NBITS { + encstate = (encstate << 1) | input[k] as u32; + coded[2 * k] = ((encstate & POLY1).count_ones() & 1) as u8; + coded[2 * k + 1] = ((encstate & POLY2).count_ones() & 1) as u8; + } + coded + } + + /// Interleave coded bits (inverse of deinterleave). + fn interleave(coded: &[u8; NSYMS]) -> [u8; NSYMS] { + let mut out = [0u8; NSYMS]; + let mut p = 0usize; + for i in 0u16..=255 { + let j = rev8(i as u8) as usize; + if j < NSYMS { + out[j] = coded[p]; + p += 1; + } + } + out + } + + /// End-to-end test: encode K1JT FN20 37, produce perfect soft symbols, + /// and verify round-trip decode. + #[test] + fn roundtrip_encode_decode() { + let cs37 = b" 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + let cs27 = b" ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + let idx37 = |c: u8| cs37.iter().position(|&x| x == c).unwrap() as u32; + let idx27 = |c: u8| cs27.iter().position(|&x| x == c).unwrap() as u32; + + let c0 = idx37(b' '); + let c1 = idx37(b'K'); + let c2 = 1u32; + let c3 = idx27(b'J'); + let c4 = idx27(b'T'); + let c5 = idx27(b' '); + let n1 = ((c0 * 36 + c1) * 10 + c2) * 27u32.pow(3) + c3 * 27u32.pow(2) + c4 * 27 + c5; + let m1 = (179 - 10 * 5 - 2) * 180 + 10 * 13 + 0; // FN20 + let power_code = 37u32; + + let mut input_bits = [0u8; NBITS]; + for i in (0..28).rev() { + input_bits[27 - i] = ((n1 >> i) & 1) as u8; + } + for i in (0..15).rev() { + input_bits[42 - i] = ((m1 >> i) & 1) as u8; + } + for i in (0..7).rev() { + input_bits[49 - i] = ((power_code >> i) & 1) as u8; + } + // bits 50..80 are tail (zeros), already set + + // Convolutional encode + let coded = convolutional_encode(&input_bits); + + // Interleave + let interleaved = interleave(&coded); + + // Create channel symbols: symbol[i] = sync[i] + 2*data_bit[i] + let channel_syms: Vec = (0..NSYMS) + .map(|i| SYNC_VECTOR[i] + 2 * interleaved[i]) + .collect(); + + // Create perfect soft symbols from channel symbols. + // data_bit = channel_sym >> 1. Soft: 0 if data=0, 255 if data=1. + let soft: Vec = channel_syms + .iter() + .map(|&cs| if cs >> 1 == 1 { 255u8 } else { 0u8 }) + .collect(); + + // Decode + let result = decode_symbols(&soft); + assert!(result.is_some(), "decode_symbols should succeed"); + let msg = result.unwrap().message; + assert!(msg.contains("K1JT"), "callsign not found in '{msg}'"); + assert!(msg.contains("FN20"), "grid not found in '{msg}'"); + assert!(msg.contains("37"), "power not found in '{msg}'"); + } + + /// Verify deinterleave is the inverse of interleave. + #[test] + fn interleave_deinterleave_roundtrip() { + // Create a sequence of distinguishable values + let mut original = [0u8; NSYMS]; + for i in 0..NSYMS { + original[i] = (i % 256) as u8; + } + + let interleaved = interleave(&original); + let recovered = deinterleave(&interleaved); + assert_eq!(original, recovered, "deinterleave should invert interleave"); + } + + /// Verify that the Fano decoder can decode a convolutionally-encoded message + /// with perfect soft symbols (0 and 255). + #[test] + fn fano_decode_perfect_soft_symbols() { + // Create a simple 81-bit message (50 payload + 31 tail zeros) + let mut input_bits = [0u8; NBITS]; + // Set some payload bits to a recognizable pattern + input_bits[0] = 1; + input_bits[5] = 1; + input_bits[10] = 1; + input_bits[15] = 1; + input_bits[20] = 1; + + // Encode + let coded = convolutional_encode(&input_bits); + + // Convert to perfect soft symbols: coded_bit=0 → 0, coded_bit=1 → 255 + let mut soft = [0u8; NSYMS]; + for i in 0..NSYMS { + soft[i] = if coded[i] == 1 { 255 } else { 0 }; + } + + // Fano decode (already in coded order, no interleaving needed) + let decoded = fano_decode(&soft); + assert!(decoded.is_some(), "Fano decoder should succeed"); + let decoded = decoded.unwrap(); + assert_eq!( + &decoded[..NBITS], + &input_bits[..NBITS], + "Decoded bits should match input" + ); + } }