mfsk_core/core/
llr.rs

1//! Protocol-agnostic soft-decision LLR computation.
2//!
3//! Extracts complex tone spectra for each data symbol, then computes four
4//! log-likelihood ratio variants (llra/b/c/d) matching WSJT-X `ft8b.f90`
5//! convention — three `nsym = 1, 2, 3` grouping schemes plus a bit-by-bit
6//! normalised variant. Parameterised over any [`Protocol`]: NTONES,
7//! BITS_PER_SYMBOL, and the SYNC_BLOCKS layout drive the inner loops.
8
9use super::Protocol;
10use num_complex::Complex;
11use rustfft::FftPlanner;
12
13// ──────────────────────────────────────────────────────────────────────────
14// LLR bundle
15// ──────────────────────────────────────────────────────────────────────────
16
17/// Four LLR vectors (a, b, c, d) of length `codeword_len()` bits in LDPC
18/// bit-index order.
19#[derive(Clone)]
20pub struct LlrSet {
21    /// nsym=1 soft metrics, scaled (matches WSJT-X llra).
22    pub llra: Vec<f32>,
23    /// nsym=2 soft metrics, scaled (matches WSJT-X llrb).
24    pub llrb: Vec<f32>,
25    /// nsym=3 soft metrics, scaled (matches WSJT-X llrc).
26    pub llrc: Vec<f32>,
27    /// nsym=1 bit-normalised (matches WSJT-X llrd).
28    pub llrd: Vec<f32>,
29}
30
31/// Default LLR scale factor from WSJT-X ft8b.f90. Individual protocols may
32/// override via `ModulationParams::LLR_SCALE`.
33pub const LLR_SCALE: f32 = 2.83;
34
35// ──────────────────────────────────────────────────────────────────────────
36// Symbol spectra
37// ──────────────────────────────────────────────────────────────────────────
38
39/// Extract complex tone spectra for every channel symbol.
40///
41/// Returns a flat row-major `Vec<Complex<f32>>` of length `N_SYMBOLS × NTONES`;
42/// row `k` / column `t` holds the k-th symbol's t-th tone amplitude, scaled
43/// by 1/1000 (matching WSJT-X).
44///
45/// `i_start` is the sample index in `cd0` of the first symbol, from fine sync.
46pub fn symbol_spectra<P: Protocol>(cd0: &[Complex<f32>], i_start: usize) -> Vec<Complex<f32>> {
47    let ntones = P::NTONES as usize;
48    let n_sym = P::N_SYMBOLS as usize;
49    let ds_spb = (P::NSPS / P::NDOWN) as usize;
50
51    let mut planner = FftPlanner::<f32>::new();
52    let fft = planner.plan_fft_forward(ds_spb);
53
54    let mut cs = vec![Complex::new(0.0f32, 0.0); n_sym * ntones];
55    let mut buf = vec![Complex::new(0.0f32, 0.0); ds_spb];
56
57    for k in 0..n_sym {
58        let i1 = i_start + k * ds_spb;
59        for (j, b) in buf.iter_mut().enumerate() {
60            *b = if i1 + j < cd0.len() {
61                cd0[i1 + j]
62            } else {
63                Complex::new(0.0, 0.0)
64            };
65        }
66        fft.process(&mut buf);
67        for (t, bin) in buf.iter().take(ntones).enumerate() {
68            cs[k * ntones + t] = *bin / 1000.0;
69        }
70    }
71    cs
72}
73
74// ──────────────────────────────────────────────────────────────────────────
75// Helpers
76// ──────────────────────────────────────────────────────────────────────────
77
78// Data-chunk layout (slots between / around sync blocks) is shared
79// with the TX side; reuse [`crate::core::tx::data_chunks`] so any
80// frame layout the encoder honours is decoded the same way.
81use crate::core::tx::data_chunks;
82
83/// Decompose `i` into `nsym` base-`ntones` digits, most significant first.
84#[inline]
85fn base_digits(mut i: usize, ntones: usize, nsym: usize) -> Vec<usize> {
86    let mut out = vec![0usize; nsym];
87    for j in (0..nsym).rev() {
88        out[j] = i % ntones;
89        i /= ntones;
90    }
91    out
92}
93
94#[inline]
95fn normalize_bmet(bmet: &mut [f32]) {
96    let n = bmet.len() as f32;
97    let mean = bmet.iter().sum::<f32>() / n;
98    let mean_sq = bmet.iter().map(|&x| x * x).sum::<f32>() / n;
99    let var = mean_sq - mean * mean;
100    let sig = if var > 0.0 {
101        var.sqrt()
102    } else {
103        mean_sq.sqrt()
104    };
105    if sig > 0.0 {
106        bmet.iter_mut().for_each(|x| *x /= sig);
107    }
108}
109
110// ──────────────────────────────────────────────────────────────────────────
111// LLR computation
112// ──────────────────────────────────────────────────────────────────────────
113
114/// Compute soft LLRs from the flat symbol-spectra vector.
115pub fn compute_llr<P: Protocol>(cs: &[Complex<f32>]) -> LlrSet {
116    let ntones = P::NTONES as usize;
117    let bps = P::BITS_PER_SYMBOL as usize;
118    let gray_map = P::GRAY_MAP;
119    let chunks = data_chunks::<P>();
120    let codeword_len: usize = chunks.iter().map(|&(_, l)| l).sum::<usize>() * bps;
121
122    let mut bmeta = vec![0.0f32; codeword_len];
123    let mut bmetb = vec![0.0f32; codeword_len];
124    let mut bmetc = vec![0.0f32; codeword_len];
125    let mut bmetd = vec![0.0f32; codeword_len];
126
127    for nsym in 1usize..=3 {
128        // Number of tone-combinations over `nsym` symbols.
129        let nt = ntones.pow(nsym as u32);
130        let ibmax = bps * nsym - 1;
131        let mut s2 = vec![0.0f32; nt];
132
133        // Walk the data symbols chunk by chunk to build the contiguous
134        // 174-bit codeword layout used by the LDPC decoder.
135        let mut chunk_bit_base = 0usize;
136        for &(chunk_start_sym, chunk_len) in &chunks {
137            let mut k = 0usize; // symbol index within this chunk
138            while k + nsym <= chunk_len {
139                let ks = chunk_start_sym + k;
140
141                // Precompute |Σ cs_k[gray[idx_k]]| for each tone-combination.
142                for (i, s2_i) in s2.iter_mut().enumerate() {
143                    let digits = base_digits(i, ntones, nsym);
144                    let sum: Complex<f32> = (0..nsym)
145                        .map(|j| cs[(ks + j) * ntones + gray_map[digits[j]] as usize])
146                        .sum();
147                    *s2_i = sum.norm();
148                }
149
150                // Map each of the `ibmax+1` bits into the codeword.
151                let i_bit_base = chunk_bit_base + k * bps;
152                for ib in 0..=ibmax {
153                    let bit_idx = i_bit_base + ib;
154                    if bit_idx >= codeword_len {
155                        break;
156                    }
157                    let bit_sel = ibmax - ib;
158                    let max_one = s2
159                        .iter()
160                        .enumerate()
161                        .filter(|&(i, _)| (i >> bit_sel) & 1 == 1)
162                        .map(|(_, &v)| v)
163                        .fold(f32::NEG_INFINITY, f32::max);
164                    let max_zero = s2
165                        .iter()
166                        .enumerate()
167                        .filter(|&(i, _)| (i >> bit_sel) & 1 == 0)
168                        .map(|(_, &v)| v)
169                        .fold(f32::NEG_INFINITY, f32::max);
170                    let bm = max_one - max_zero;
171
172                    match nsym {
173                        1 => {
174                            bmeta[bit_idx] = bm;
175                            let den = max_one.max(max_zero);
176                            bmetd[bit_idx] = if den > 0.0 { bm / den } else { 0.0 };
177                        }
178                        2 => bmetb[bit_idx] = bm,
179                        3 => bmetc[bit_idx] = bm,
180                        _ => unreachable!(),
181                    }
182                }
183
184                k += nsym;
185            }
186            chunk_bit_base += chunk_len * bps;
187        }
188    }
189
190    normalize_bmet(&mut bmeta);
191    normalize_bmet(&mut bmetb);
192    normalize_bmet(&mut bmetc);
193    normalize_bmet(&mut bmetd);
194
195    let s = P::LLR_SCALE;
196    let scale = |v: Vec<f32>| -> Vec<f32> { v.into_iter().map(|x| x * s).collect() };
197
198    LlrSet {
199        llra: scale(bmeta),
200        llrb: scale(bmetb),
201        llrc: scale(bmetc),
202        llrd: scale(bmetd),
203    }
204}
205
206// ──────────────────────────────────────────────────────────────────────────
207// SNR estimation
208// ──────────────────────────────────────────────────────────────────────────
209
210/// WSJT-X compatible SNR (dB) estimate from symbol spectra + decoded tones.
211///
212/// Signal: `Σ |cs[k][itone[k]]|²`. Noise reference: `Σ |cs[k][(itone[k] +
213/// NTONES/2) mod NTONES]|²` (tone on the "opposite side" of the comb).
214/// SNR_dB = `10·log10(sig/noi − 1) − 27` clamped to −24 dB floor (WSJT-X
215/// convention, applied per-tone bandwidth → 2500 Hz reference).
216pub fn compute_snr_db<P: Protocol>(cs: &[Complex<f32>], itone: &[u8]) -> f32 {
217    let ntones = P::NTONES as usize;
218    let n_sym = P::N_SYMBOLS as usize;
219    let mut xsig = 0.0f32;
220    let mut xnoi = 0.0f32;
221    let offset = ntones / 2;
222    for k in 0..n_sym.min(itone.len()) {
223        let t = itone[k] as usize % ntones;
224        xsig += cs[k * ntones + t].norm_sqr();
225        xnoi += cs[k * ntones + (t + offset) % ntones].norm_sqr();
226    }
227    if xnoi < f32::EPSILON {
228        return -24.0;
229    }
230    let ratio = xsig / xnoi - 1.0;
231    if ratio <= 0.001 {
232        return -24.0;
233    }
234    (10.0 * ratio.log10() - 27.0_f32).max(-24.0)
235}
236
237/// Hard-decision sync quality — count sync symbols whose dominant tone
238/// matches the protocol's Costas pattern. Range is 0..N_SYNC; callers
239/// typically threshold on this.
240pub fn sync_quality<P: Protocol>(cs: &[Complex<f32>]) -> u32 {
241    let ntones = P::NTONES as usize;
242    let mut count = 0u32;
243    for block in P::SYNC_MODE.blocks() {
244        let start = block.start_symbol as usize;
245        for (t, &expected) in block.pattern.iter().enumerate() {
246            let sym = start + t;
247            let best = (0..ntones)
248                .max_by(|&a, &b| {
249                    cs[sym * ntones + a]
250                        .norm()
251                        .partial_cmp(&cs[sym * ntones + b].norm())
252                        .unwrap()
253                })
254                .unwrap_or(0);
255            if best == expected as usize {
256                count += 1;
257            }
258        }
259    }
260    count
261}
mfsk_core/core/llr.rs

mfsk_core/core/
llr.rs
