Signed-off-by: Stan Grams <sjg@haxx.space>
This commit is contained in:
@@ -0,0 +1,10 @@
|
||||
# SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
#
|
||||
# SPDX-License-Identifier: GPL-2.0-or-later
|
||||
|
||||
[package]
|
||||
name = "trx-wspr"
|
||||
version.workspace = true
|
||||
edition = "2021"
|
||||
|
||||
[dependencies]
|
||||
@@ -0,0 +1,510 @@
|
||||
// SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
//
|
||||
// SPDX-License-Identifier: GPL-2.0-or-later
|
||||
|
||||
use crate::protocol;
|
||||
|
||||
const WSPR_SAMPLE_RATE: u32 = 12_000;
|
||||
const SLOT_SAMPLES: usize = 120 * WSPR_SAMPLE_RATE as usize;
|
||||
const WSPR_SYMBOL_COUNT: usize = 162;
|
||||
const WSPR_SYMBOL_SAMPLES: usize = 8192;
|
||||
const WSPR_SIGNAL_SAMPLES: usize = WSPR_SYMBOL_COUNT * WSPR_SYMBOL_SAMPLES;
|
||||
const EXPECTED_SIGNAL_START_SAMPLES: usize = WSPR_SAMPLE_RATE as usize; // 1s
|
||||
const TONE_SPACING_HZ: f32 = WSPR_SAMPLE_RATE as f32 / WSPR_SYMBOL_SAMPLES as f32; // 1.46484375
|
||||
|
||||
// Coarse search range for base tone. This matches common WSPR audio passband.
|
||||
const BASE_SEARCH_MIN_HZ: f32 = 1200.0;
|
||||
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 (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 normalized sync correlation score to attempt decode.
|
||||
// The reference wsprd uses minsync1=0.10 but applies additional filtering
|
||||
// downstream. A higher threshold here prevents noise from reaching the Fano
|
||||
// decoder and producing false positives.
|
||||
const MIN_SYNC_SCORE: f32 = 0.20;
|
||||
|
||||
// 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]
|
||||
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,
|
||||
1,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1,
|
||||
0,1,0,0,0,1,1,1,0,0,0,0,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,1,1,0,
|
||||
1,0,1,1,0,0,0,1,1,0,0,1,
|
||||
];
|
||||
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct WsprDecodeResult {
|
||||
pub message: String,
|
||||
pub snr_db: f32,
|
||||
pub dt_s: f32,
|
||||
pub freq_hz: f32,
|
||||
}
|
||||
|
||||
// Minimum estimated SNR (dB) to attempt decode. WSPR's theoretical decode
|
||||
// limit is around -28 dB in 2500 Hz bandwidth, but the per-tone SNR estimate
|
||||
// computed here uses a narrower noise reference and reads higher. Setting
|
||||
// -20 dB is conservative enough to pass all real signals while rejecting
|
||||
// pure-noise candidates where the Fano decoder might otherwise hallucinate.
|
||||
const MIN_SNR_DB: f32 = -20.0;
|
||||
|
||||
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.005 })
|
||||
}
|
||||
|
||||
pub fn sample_rate(&self) -> u32 {
|
||||
WSPR_SAMPLE_RATE
|
||||
}
|
||||
|
||||
pub fn slot_samples(&self) -> usize {
|
||||
SLOT_SAMPLES
|
||||
}
|
||||
|
||||
pub fn decode_slot(
|
||||
&self,
|
||||
samples: &[f32],
|
||||
_base_freq_hz: Option<u64>,
|
||||
) -> Result<Vec<WsprDecodeResult>, String> {
|
||||
if samples.len() < SLOT_SAMPLES {
|
||||
return Ok(Vec::new());
|
||||
}
|
||||
|
||||
let rms = slot_rms(&samples[..SLOT_SAMPLES]);
|
||||
if rms < self.min_rms {
|
||||
return Ok(Vec::new());
|
||||
}
|
||||
|
||||
// Collect top frequency candidates across timing offsets
|
||||
let mut candidates: Vec<(f32, isize, f32)> = Vec::new(); // (freq, dt_samples, score)
|
||||
|
||||
let mut dt = -DT_SEARCH_RANGE_SAMPLES;
|
||||
while dt <= DT_SEARCH_RANGE_SAMPLES {
|
||||
let start = EXPECTED_SIGNAL_START_SAMPLES as isize + dt;
|
||||
if start < 0 || (start as usize) + WSPR_SIGNAL_SAMPLES > samples.len() {
|
||||
dt += DT_SEARCH_STEP_SAMPLES;
|
||||
continue;
|
||||
}
|
||||
let signal = &samples[start as usize..start as usize + WSPR_SIGNAL_SAMPLES];
|
||||
|
||||
// Coarse frequency search using sync vector correlation
|
||||
let mut freq_scores: Vec<(f32, f32)> = Vec::new();
|
||||
let mut freq = BASE_SEARCH_MIN_HZ;
|
||||
while freq <= BASE_SEARCH_MAX_HZ {
|
||||
let score = sync_correlation_score(signal, freq);
|
||||
freq_scores.push((freq, score));
|
||||
freq += BASE_SEARCH_STEP_HZ;
|
||||
}
|
||||
|
||||
// Keep top candidates from coarse search
|
||||
freq_scores.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
|
||||
for &(coarse_freq, _) in freq_scores.iter().take(3) {
|
||||
// Fine-tune frequency around each coarse candidate
|
||||
let mut best_fine_freq = coarse_freq;
|
||||
let mut best_fine_score = f32::MIN;
|
||||
let mut fine_freq = coarse_freq - BASE_SEARCH_STEP_HZ;
|
||||
while fine_freq <= coarse_freq + BASE_SEARCH_STEP_HZ {
|
||||
let score = sync_correlation_score(signal, fine_freq);
|
||||
if score > best_fine_score {
|
||||
best_fine_score = score;
|
||||
best_fine_freq = fine_freq;
|
||||
}
|
||||
fine_freq += FINE_SEARCH_STEP_HZ;
|
||||
}
|
||||
candidates.push((best_fine_freq, dt, best_fine_score));
|
||||
}
|
||||
dt += DT_SEARCH_STEP_SAMPLES;
|
||||
}
|
||||
|
||||
// Sort candidates by score (best first) and try to decode each
|
||||
candidates.sort_by(|a, b| b.2.partial_cmp(&a.2).unwrap_or(std::cmp::Ordering::Equal));
|
||||
|
||||
let mut results = Vec::new();
|
||||
let mut seen_messages = std::collections::HashSet::new();
|
||||
// Track (freq, dt) of successful decodes to skip near-duplicates
|
||||
let mut decoded_positions: Vec<(f32, isize)> = Vec::new();
|
||||
|
||||
for &(freq, dt_samples, _score) in candidates.iter().take(MAX_FREQ_CANDIDATES) {
|
||||
// Skip candidates too close in (freq, dt) to an already-decoded signal
|
||||
let dominated = decoded_positions.iter().any(|&(df, ddt)| {
|
||||
(freq - df).abs() < 4.0 * TONE_SPACING_HZ
|
||||
&& (dt_samples - ddt).unsigned_abs() < DT_SEARCH_STEP_SAMPLES as usize
|
||||
});
|
||||
if dominated {
|
||||
continue;
|
||||
}
|
||||
let start = (EXPECTED_SIGNAL_START_SAMPLES as isize + dt_samples) as usize;
|
||||
let signal = &samples[start..start + WSPR_SIGNAL_SAMPLES];
|
||||
|
||||
// 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);
|
||||
|
||||
// Reject candidates where estimated SNR is too low — the Fano
|
||||
// decoder can converge on noise-only input after normalization.
|
||||
if demod.snr_db < MIN_SNR_DB {
|
||||
continue;
|
||||
}
|
||||
|
||||
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 {
|
||||
message: decoded.message,
|
||||
snr_db: demod.snr_db,
|
||||
dt_s,
|
||||
freq_hz: freq,
|
||||
});
|
||||
decoded_positions.push((freq, dt_samples));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Ok(results)
|
||||
}
|
||||
}
|
||||
|
||||
/// 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 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];
|
||||
|
||||
// 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: (p1+p3)-(p0+p2) weighted by (2*sync-1)
|
||||
let cmet = (p1 + p3) - (p0 + p2);
|
||||
if sync_bit == 1 {
|
||||
ss += cmet;
|
||||
} else {
|
||||
ss -= cmet;
|
||||
}
|
||||
}
|
||||
|
||||
// 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
|
||||
}
|
||||
|
||||
/// 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;
|
||||
|
||||
for sym in 0..WSPR_SYMBOL_COUNT {
|
||||
let off = sym * WSPR_SYMBOL_SAMPLES;
|
||||
let frame = &signal[off..off + WSPR_SYMBOL_SAMPLES];
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
// 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;
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
// 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::<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 {
|
||||
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();
|
||||
let w = 2.0 * std::f32::consts::PI * k / n;
|
||||
let coeff = 2.0 * w.cos();
|
||||
|
||||
let mut s_prev = 0.0_f32;
|
||||
let mut s_prev2 = 0.0_f32;
|
||||
for (idx, &x) in frame.iter().enumerate() {
|
||||
let win = 0.5_f32 - 0.5_f32 * (2.0_f32 * std::f32::consts::PI * idx as f32 / n).cos();
|
||||
let s = x * win + coeff * s_prev - s_prev2;
|
||||
s_prev2 = s_prev;
|
||||
s_prev = s;
|
||||
}
|
||||
|
||||
s_prev2 * s_prev2 + s_prev * s_prev - coeff * s_prev * s_prev2
|
||||
}
|
||||
|
||||
fn slot_rms(samples: &[f32]) -> f32 {
|
||||
if samples.is_empty() {
|
||||
return 0.0;
|
||||
}
|
||||
let sum_sq = samples.iter().map(|s| s * s).sum::<f32>();
|
||||
(sum_sq / samples.len() as f32).sqrt()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn short_slot_returns_empty() {
|
||||
let dec = WsprDecoder::new().expect("decoder");
|
||||
let out = dec.decode_slot(&vec![0.0; dec.slot_samples() - 1], None);
|
||||
assert!(out.expect("decode").is_empty());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn rms_is_zero_for_silence() {
|
||||
let rms = slot_rms(&[0.0; 16]);
|
||||
assert_eq!(rms, 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn base_search_finds_synthetic_signal() {
|
||||
let mut slot = vec![0.0_f32; SLOT_SAMPLES];
|
||||
let base_hz = 1496.0_f32;
|
||||
let start = EXPECTED_SIGNAL_START_SAMPLES;
|
||||
|
||||
for (sym, sync_tone) in SYNC_VECTOR
|
||||
.iter()
|
||||
.copied()
|
||||
.enumerate()
|
||||
.take(WSPR_SYMBOL_COUNT)
|
||||
{
|
||||
let tone = sync_tone + 2 * ((sym % 2) as u8);
|
||||
let freq = base_hz + tone as f32 * TONE_SPACING_HZ;
|
||||
let begin = start + sym * WSPR_SYMBOL_SAMPLES;
|
||||
for i in 0..WSPR_SYMBOL_SAMPLES {
|
||||
let t = i as f32 / WSPR_SAMPLE_RATE as f32;
|
||||
slot[begin + i] = (2.0 * std::f32::consts::PI * freq * t).sin() * 0.2;
|
||||
}
|
||||
}
|
||||
|
||||
let signal = &slot[start..start + WSPR_SIGNAL_SAMPLES];
|
||||
let candidates = find_candidates(signal);
|
||||
assert!(!candidates.is_empty());
|
||||
let (estimated, _) = candidates[0];
|
||||
assert!(
|
||||
(estimated - base_hz).abs() <= 1.0,
|
||||
"estimated {estimated} Hz, expected {base_hz} Hz"
|
||||
);
|
||||
}
|
||||
|
||||
/// Helper: run the candidate search on a signal slice
|
||||
fn find_candidates(signal: &[f32]) -> Vec<(f32, f32)> {
|
||||
let mut freq_scores: Vec<(f32, f32)> = Vec::new();
|
||||
let mut freq = BASE_SEARCH_MIN_HZ;
|
||||
while freq <= BASE_SEARCH_MAX_HZ {
|
||||
let score = sync_correlation_score(signal, freq);
|
||||
freq_scores.push((freq, score));
|
||||
freq += BASE_SEARCH_STEP_HZ;
|
||||
}
|
||||
freq_scores.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
|
||||
|
||||
// Fine-tune top result
|
||||
if let Some(&(coarse_freq, _)) = freq_scores.first() {
|
||||
let mut best_fine_freq = coarse_freq;
|
||||
let mut best_fine_score = f32::MIN;
|
||||
let mut fine_freq = coarse_freq - BASE_SEARCH_STEP_HZ;
|
||||
while fine_freq <= coarse_freq + BASE_SEARCH_STEP_HZ {
|
||||
let score = sync_correlation_score(signal, fine_freq);
|
||||
if score > best_fine_score {
|
||||
best_fine_score = score;
|
||||
best_fine_freq = fine_freq;
|
||||
}
|
||||
fine_freq += FINE_SEARCH_STEP_HZ;
|
||||
}
|
||||
vec![(best_fine_freq, best_fine_score)]
|
||||
} else {
|
||||
vec![]
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sync_correlation_prefers_correct_frequency() {
|
||||
let base_hz = 1500.0_f32;
|
||||
let wrong_hz = 1400.0_f32;
|
||||
|
||||
// Generate a synthetic WSPR-like signal using the sync vector
|
||||
let mut signal = vec![0.0_f32; WSPR_SIGNAL_SAMPLES];
|
||||
for (sym, sync_tone) in SYNC_VECTOR
|
||||
.iter()
|
||||
.copied()
|
||||
.enumerate()
|
||||
.take(WSPR_SYMBOL_COUNT)
|
||||
{
|
||||
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 correct_score = sync_correlation_score(&signal, base_hz);
|
||||
let wrong_score = sync_correlation_score(&signal, wrong_hz);
|
||||
assert!(
|
||||
correct_score > wrong_score,
|
||||
"correct={correct_score}, wrong={wrong_score}"
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noise_only_slot_produces_no_decodes() {
|
||||
// Deterministic pseudo-random noise via simple LCG
|
||||
let mut rng_state = 0x12345678u64;
|
||||
let mut next_f32 = || -> f32 {
|
||||
rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
|
||||
((rng_state >> 33) as f32 / u32::MAX as f32) * 2.0 - 1.0
|
||||
};
|
||||
|
||||
let dec = WsprDecoder::new().expect("decoder");
|
||||
let slot: Vec<f32> = (0..dec.slot_samples()).map(|_| next_f32() * 0.05).collect();
|
||||
let results = dec.decode_slot(&slot, None).expect("decode");
|
||||
assert!(
|
||||
results.is_empty(),
|
||||
"noise-only slot should produce no decodes, got {}",
|
||||
results.len()
|
||||
);
|
||||
}
|
||||
|
||||
#[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}");
|
||||
// This synthetic signal only uses sync tones (no data tones), so the
|
||||
// normalized score is moderate (~0.18). A real WSPR signal occupies all
|
||||
// 4 tones and produces higher scores (>0.3).
|
||||
assert!(
|
||||
score > 0.10,
|
||||
"score {score} should be clearly above noise floor"
|
||||
);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,8 @@
|
||||
// SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
//
|
||||
// SPDX-License-Identifier: GPL-2.0-or-later
|
||||
|
||||
mod decoder;
|
||||
mod protocol;
|
||||
|
||||
pub use decoder::{WsprDecodeResult, WsprDecoder};
|
||||
@@ -0,0 +1,578 @@
|
||||
// SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
//
|
||||
// SPDX-License-Identifier: GPL-2.0-or-later
|
||||
|
||||
/// Decoded WSPR message payload.
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct WsprProtocolMessage {
|
||||
pub message: String,
|
||||
}
|
||||
|
||||
const POLY1: u32 = 0xF2D05351;
|
||||
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;
|
||||
for _ in 0..8 {
|
||||
r = (r << 1) | (b & 1);
|
||||
b >>= 1;
|
||||
}
|
||||
r
|
||||
}
|
||||
|
||||
/// Deinterleave soft symbols by permuting their order via bit-reversal of indices.
|
||||
///
|
||||
/// 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 = [128u8; NSYMS]; // default to "no confidence"
|
||||
let mut p = 0usize;
|
||||
for i in 0u16..=255 {
|
||||
let j = rev8(i as u8) as usize;
|
||||
if j < NSYMS {
|
||||
out[p] = if j < symbols.len() { symbols[j] } else { 128 };
|
||||
p += 1;
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// Compute the 2-bit convolutional encoder output for a given encoder state.
|
||||
///
|
||||
/// 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
|
||||
}
|
||||
|
||||
/// Result from the Fano decoder including quality metric.
|
||||
struct FanoResult {
|
||||
bits: [u8; NBITS],
|
||||
/// Cumulative path metric — higher values indicate higher confidence.
|
||||
metric: i64,
|
||||
}
|
||||
|
||||
/// 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: decoded bits and cumulative path metric, or None on timeout.
|
||||
fn fano_decode(symbols: &[u8; NSYMS]) -> Option<FanoResult> {
|
||||
let mettab = build_mettab();
|
||||
let max_cycles = FANO_MAX_CYCLES_PER_BIT * NBITS;
|
||||
let tail_start = NBITS - 31; // position 50: first tail bit
|
||||
|
||||
// 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 t: i64 = 0; // threshold
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
// 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; // 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(FanoResult {
|
||||
bits,
|
||||
metric: gamma[NBITS],
|
||||
})
|
||||
}
|
||||
|
||||
/// Unpack 50 payload bits into a formatted WSPR message string.
|
||||
///
|
||||
/// Layout (MSB first):
|
||||
/// bits 0-27 — N1 (28 bits): callsign
|
||||
/// bits 28-42 — M1 (15 bits): Maidenhead grid
|
||||
/// bits 43-49 — P ( 7 bits): power code (dBm + 64)
|
||||
fn unpack_message(bits: &[u8; NBITS]) -> Option<String> {
|
||||
// Accumulate N1, M1, and power code from the bit array.
|
||||
let mut n1 = 0u32;
|
||||
for &b in &bits[..28] {
|
||||
n1 = (n1 << 1) | b as u32;
|
||||
}
|
||||
let mut m1 = 0u32;
|
||||
for &b in &bits[28..43] {
|
||||
m1 = (m1 << 1) | b as u32;
|
||||
}
|
||||
let mut power_code = 0u32;
|
||||
for &b in &bits[43..50] {
|
||||
power_code = (power_code << 1) | b as u32;
|
||||
}
|
||||
|
||||
// WSPR only permits specific power levels (dBm).
|
||||
const VALID_POWER: [i32; 19] = [
|
||||
0, 3, 7, 10, 13, 17, 20, 23, 27, 30, 33, 37, 40, 43, 47, 50, 53, 57, 60,
|
||||
];
|
||||
let power_dbm = power_code as i32;
|
||||
if !VALID_POWER.contains(&power_dbm) {
|
||||
return None;
|
||||
}
|
||||
|
||||
// Decode callsign from N1.
|
||||
// N1 = ((c0*36 + c1)*10 + c2)*27^3 + c3*27^2 + c4*27 + c5
|
||||
// c0,c1 ∈ charset37; c2 ∈ '0'-'9'; c3,c4,c5 ∈ charset27
|
||||
const CS37: &[u8] = b" 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
|
||||
const CS27: &[u8] = b" ABCDEFGHIJKLMNOPQRSTUVWXYZ";
|
||||
|
||||
let mut n = n1;
|
||||
let i5 = (n % 27) as usize;
|
||||
n /= 27;
|
||||
let i4 = (n % 27) as usize;
|
||||
n /= 27;
|
||||
let i3 = (n % 27) as usize;
|
||||
n /= 27;
|
||||
let i2 = (n % 10) as usize;
|
||||
n /= 10;
|
||||
let i1 = (n % 36) as usize;
|
||||
n /= 36;
|
||||
let i0 = n as usize;
|
||||
|
||||
if i0 >= 37 || i1 >= 37 || i2 >= 10 || i3 >= 27 || i4 >= 27 || i5 >= 27 {
|
||||
return None;
|
||||
}
|
||||
|
||||
let callsign = format!(
|
||||
"{}{}{}{}{}{}",
|
||||
CS37[i0] as char,
|
||||
CS37[i1] as char,
|
||||
(b'0' + i2 as u8) as char,
|
||||
CS27[i3] as char,
|
||||
CS27[i4] as char,
|
||||
CS27[i5] as char,
|
||||
)
|
||||
.trim()
|
||||
.to_string();
|
||||
|
||||
// WSPR callsigns: after trimming, the digit (from position 2 of the
|
||||
// 6-char padded form) must appear at index 1 or 2. The callsign must
|
||||
// also contain at least one letter and be at least 3 characters long.
|
||||
if callsign.len() < 3 || !callsign.chars().any(|c| c.is_alphabetic()) {
|
||||
return None;
|
||||
}
|
||||
let has_digit_at_1_or_2 = callsign.chars().nth(1).is_some_and(|c| c.is_ascii_digit())
|
||||
|| callsign.chars().nth(2).is_some_and(|c| c.is_ascii_digit());
|
||||
if !has_digit_at_1_or_2 {
|
||||
return None;
|
||||
}
|
||||
|
||||
// Decode Maidenhead grid from M1.
|
||||
// M1 = (179 - 10*loc1 - loc3)*180 + 10*loc2 + loc4
|
||||
// loc1,loc2 ∈ 0-17 (A-R); loc3,loc4 ∈ 0-9
|
||||
if m1 > 32_399 {
|
||||
return None;
|
||||
}
|
||||
let hi = m1 / 180;
|
||||
let lo = m1 % 180;
|
||||
let t = 179u32.checked_sub(hi)?;
|
||||
let loc1 = t / 10; // longitude letter index
|
||||
let loc3 = t % 10; // longitude digit
|
||||
let loc2 = lo / 10; // latitude letter index
|
||||
let loc4 = lo % 10; // latitude digit
|
||||
|
||||
if loc1 > 17 || loc2 > 17 {
|
||||
return None;
|
||||
}
|
||||
|
||||
let grid = format!(
|
||||
"{}{}{}{}",
|
||||
(b'A' + loc1 as u8) as char,
|
||||
(b'A' + loc2 as u8) as char,
|
||||
(b'0' + loc3 as u8) as char,
|
||||
(b'0' + loc4 as u8) as char,
|
||||
);
|
||||
|
||||
Some(format!("{} {} {}", callsign, grid, power_dbm))
|
||||
}
|
||||
|
||||
/// Minimum Fano cumulative path metric to accept a decode.
|
||||
///
|
||||
/// The Fano decoder can sometimes converge on random noise, producing bits
|
||||
/// that happen to unpack into a valid-looking message. The cumulative path
|
||||
/// metric reflects how well the received symbols matched the best trellis
|
||||
/// path. Real WSPR signals at decodable SNR produce metrics well above this
|
||||
/// threshold; noise-induced decodes have metrics near or below zero.
|
||||
const FANO_MIN_METRIC: i64 = 20;
|
||||
|
||||
/// 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;
|
||||
}
|
||||
let coded = deinterleave(symbols);
|
||||
let result = fano_decode(&coded)?;
|
||||
|
||||
// Reject low-confidence decodes that are likely false positives from noise
|
||||
if result.metric < FANO_MIN_METRIC {
|
||||
return None;
|
||||
}
|
||||
|
||||
let message = unpack_message(&result.bits)?;
|
||||
Some(WsprProtocolMessage { message })
|
||||
}
|
||||
|
||||
#[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.
|
||||
#[test]
|
||||
fn unpack_known_message() {
|
||||
// Callsign "K1JT", grid "FN20", power 37 dBm — a well-known WSPR beacon.
|
||||
// Encode callsign "K1JT " (padded to 6 chars with trailing spaces).
|
||||
// charset37: ' '=0, '0'=1,..'9'=10, 'A'=11,..'Z'=36
|
||||
// charset27: ' '=0, 'A'=1,..'Z'=26
|
||||
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;
|
||||
|
||||
// " 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');
|
||||
let c2 = 1u32; // '1'
|
||||
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;
|
||||
|
||||
// Grid "FN20": loc1='F'=5 (lon), loc2='N'=13 (lat), loc3='2', loc4='0'
|
||||
let loc1 = (b'F' - b'A') as u32; // 5
|
||||
let loc2 = (b'N' - b'A') as u32; // 13
|
||||
let loc3 = 2u32;
|
||||
let loc4 = 0u32;
|
||||
let m1 = (179 - 10 * loc1 - loc3) * 180 + 10 * loc2 + loc4;
|
||||
|
||||
// Power 37 dBm → power_code = 37 (raw dBm value)
|
||||
let power_code = 37u32;
|
||||
|
||||
// Pack into 50-bit array
|
||||
let mut bits = [0u8; NBITS];
|
||||
for i in (0..28).rev() {
|
||||
bits[27 - i] = ((n1 >> i) & 1) as u8;
|
||||
}
|
||||
for i in (0..15).rev() {
|
||||
bits[42 - i] = ((m1 >> i) & 1) as u8;
|
||||
}
|
||||
for i in (0..7).rev() {
|
||||
bits[49 - i] = ((power_code >> i) & 1) as u8;
|
||||
}
|
||||
|
||||
let msg = unpack_message(&bits).expect("unpack_message should succeed");
|
||||
// Message should contain callsign, grid, and power
|
||||
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);
|
||||
}
|
||||
|
||||
/// 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 result = fano_decode(&soft);
|
||||
assert!(result.is_some(), "Fano decoder should succeed");
|
||||
let result = result.unwrap();
|
||||
assert_eq!(
|
||||
&result.bits[..NBITS],
|
||||
&input_bits[..NBITS],
|
||||
"Decoded bits should match input"
|
||||
);
|
||||
assert!(
|
||||
result.metric > 0,
|
||||
"Path metric should be positive for perfect symbols"
|
||||
);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user