[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 <noreply@anthropic.com>
Signed-off-by: Stan Grams <sjg@haxx.space>

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
+// Minimum normalized sync correlation score to attempt decode.
-// this threshold are almost certainly noise and skipping them avoids expensive
+// Matches reference wsprd minsync1=0.10.
-// Fano decode attempts that would produce false positives.
+const MIN_SYNC_SCORE: f32 = 0.10;
-const MIN_SYNC_SCORE: f32 = 10.0;
+
+// 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<Self, String> {
         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) {
+        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
-            }
             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);
+            // Use normalized sync score for threshold check
-            if let Some(decoded) = protocol::decode_symbols(&demod.symbols) {
+            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)]
+/// Score a candidate base frequency by correlating detected tone amplitudes
-struct DemodOutput {
+/// with the known WSPR sync vector. Uses amplitude (sqrt of power) and
-    symbols: Vec<u8>,
+/// normalizes by total power, matching the reference wsprd implementation.
-    snr_db: f32,
+/// Higher score = better match. Range approximately [0.0, 1.0].
-}
-
-/// Score a candidate base frequency by correlating detected symbol LSBs with
-/// the known WSPR sync vector. Higher score = better match.
 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
+        // Compute amplitude (sqrt of power) at each of the 4 FSK tones
-        let odd_power = p1 + p3; // tones with LSB=1
+        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 {
+/// Compute the normalized sync score (ss/totp) for threshold comparison.
-    let mut symbols = Vec::with_capacity(WSPR_SYMBOL_COUNT);
+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];
+        // Compute amplitude (sqrt of power) at each tone — matches reference
-        for (i, power) in tone_power.iter_mut().enumerate() {
+        let p0 = goertzel_power(frame, base_hz, sr).sqrt();
-            let hz = base_hz + i as f32 * TONE_SPACING_HZ;
+        let p1 = goertzel_power(frame, base_hz + TONE_SPACING_HZ, sr).sqrt();
-            *power = goertzel_power(frame, hz, WSPR_SAMPLE_RATE as f32);
+        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;
+        // SNR estimation: signal = best tone power, noise = out-of-band
-        let mut best_pow = tone_power[0];
+        let best_amp = p0.max(p1).max(p2).max(p3);
-        for (idx, p) in tone_power.iter().enumerate().skip(1) {
+        signal_sum += best_amp * best_amp;
-            if *p > best_pow {
-                best_pow = *p;
-                best_idx = idx as u8;
-            }
-        }
 
-        symbols.push(best_idx);
+        let noise_a = goertzel_power(frame, base_hz - 8.0 * TONE_SPACING_HZ, sr);
-        signal_sum += best_pow;
+        let noise_b = goertzel_power(frame, base_hz + 12.0 * TONE_SPACING_HZ, sr);
-
-        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,
-        );
         noise_sum += (noise_a + noise_b) * 0.5;
     }
 
-    let signal_avg = signal_sum / WSPR_SYMBOL_COUNT as f32;
+    // Normalize: zero-mean, scale by symfac/stddev, clip to [-128,127], bias to [0,255]
-    let noise_avg = (noise_sum / WSPR_SYMBOL_COUNT as f32).max(1e-12);
+    let n = WSPR_SYMBOL_COUNT as f32;
+    let mean = fsymb.iter().sum::<f32>() / n;
+    let var = fsymb.iter().map(|&x| (x - mean) * (x - mean)).sum::<f32>() / 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}"
+        );
+    }
 }

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.
+/// Unlike the old hard-decision version, this does NOT extract data bits — the
-/// The 162 data bits are reordered via bit-reversal of 8-bit indices.
+/// soft values (0-255, centered at 128) are preserved as-is. The Fano decoder
-///
+/// interprets them directly via the metric table.
-/// 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.
 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 {
+    for i in 0u16..=255 {
-        let j = rev8(i) as usize;
+        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 a value 0-3 where:
-/// Returns the 81 decoded bits (first 50 are payload), or None if it cannot
+///   bit 1 (2's place) = parity(state & POLY1)
-/// converge within the iteration budget.
+///   bit 0 (1's place) = parity(state & POLY2)
-fn fano_decode(coded: &[u8; NSYMS]) -> Option<[u8; NBITS]> {
+fn encode_sym(state: u32) -> u32 {
-    const MAX_CYCLES: usize = 100_000;
+    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.
+/// Soft-decision Fano sequential decoder for K=32, rate-1/2 convolutional code.
-    let bm = |pos: usize, state: u32, bit: u32| -> i32 {
+///
-        let ns = (state << 1) | bit;
+/// Closely follows the reference implementation from WSJT-X (fano.c by Phil Karn, KA9Q).
-        let c0 = ((ns & POLY1).count_ones() & 1) as u8;
+///
-        let c1 = ((ns & POLY2).count_ones() & 1) as u8;
+/// Input: 162 deinterleaved soft-decision symbols (0-255, 128=no confidence).
-        let r0 = coded[2 * pos];
+/// Symbols are read in pairs: `symbols[2k]` and `symbols[2k+1]` are the two
-        let r1 = coded[2 * pos + 1];
+/// coded bits for input bit k.
-        (if c0 == r0 { 1 } else { -1 }) + (if c1 == r1 { 1 } else { -1 })
+///
-    };
+/// 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];
+    // Precompute all 4 branch metrics for each bit position.
-    let mut states = [0u32; NBITS + 1];
+    // metrics[k][sym_pair] where sym_pair encodes (expected_bit0, expected_bit1):
-    let mut metrics = [0i32; NBITS + 1];
+    //   0 = (0,0), 1 = (0,1), 2 = (1,0), 3 = (1,1)
-    // 0 = not yet visited, 1 = tried best branch, 2 = tried both branches
+    let mut metrics = [[0i32; 4]; NBITS];
-    let mut visit = [0u8; 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 {
+    // Initialize root node: compute and sort branch metrics
-        if pos == NBITS {
+    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];
+        // Look forward: try current branch
-        let base = metrics[pos];
+        let ngamma = gamma[pos] + tm[pos][branch_i[pos] as usize] as i64;
-        let m0 = base + bm(pos, state, 0);
+        if ngamma >= t {
-        let m1 = base + bm(pos, state, 1);
+            // Acceptable — tighten threshold if this is a first visit
-
+            if gamma[pos] < t + FANO_DELTA as i64 {
-        // b0/mv0 = best branch; b1/mv1 = second-best
+                while ngamma >= t + FANO_DELTA as i64 {
-        let (b0, mv0, b1, mv1) = if m0 >= m1 {
+                    t += FANO_DELTA as i64;
-            (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.
             }
+
+            // 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<String> {
     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<WsprProtocolMessage> {
     if symbols.len() < NSYMS {
         return None;
@@ -240,6 +362,7 @@ pub fn decode_symbols(symbols: &[u8]) -> Option<WsprProtocolMessage> {
 #[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<u8> = (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<u8> = 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"
+        );
+    }
 }