mfsk_core/core/dsp/
downsample.rs

1//! FFT-based downconversion / decimation.
2//!
3//! Ported from WSJT-X `ft8_downsample.f90` and generalised to arbitrary MFSK
4//! signals. The algorithm is:
5//!
6//! 1. Zero-pad the real input to [`DownsampleCfg::fft1_size`] samples and take
7//!    a forward FFT.
8//! 2. Extract the positive-frequency bins covering
9//!    `[f0 - leading_pad·Δf, f0 + trailing_pad·Δf]`, where the pads are
10//!    expressed in tone-spacing units.
11//! 3. Hann-taper [`DownsampleCfg::edge_taper_bins`] on each side of the
12//!    extracted block to suppress FFT ringing.
13//! 4. Cyclic-rotate so that `f0` lands at DC.
14//! 5. Inverse FFT of size [`DownsampleCfg::fft2_size`] gives the decimated
15//!    complex baseband.
16//! 6. Scale by `1 / sqrt(fft1_size · fft2_size)`.
17//!
18//! The output sample rate is `input_rate · fft2_size / fft1_size`. For the FT8
19//! tuning (12 000 → 200 Hz) it is 200 Hz; FT4 will pick a different
20//! `fft2_size` to keep a wider baseband.
21//!
22//! ## Why a config struct rather than `<P>` ?
23//!
24//! `fft1_size` is chosen per protocol for FFT efficiency (highly-composite
25//! numbers close to the slot length). It is not simply derived from
26//! `SLOT_S · sample_rate` — it is a tunable tied to the FFT backend's
27//! strengths. Keeping it in a runtime struct lets callers express values that
28//! are awkward to pin to associated constants.
29
30use alloc::vec;
31use alloc::vec::Vec;
32use core::iter;
33
34use num_complex::Complex;
35#[cfg(not(feature = "std"))]
36use num_traits::Float;
37
38use crate::core::fft::default_planner;
39
40/// Runtime parameters shared by [`downsample`], [`downsample_cached`], and
41/// [`build_fft_cache`]. Callers typically keep one instance per protocol.
42#[derive(Clone, Copy, Debug)]
43pub struct DownsampleCfg {
44    /// Input sample rate in Hz (12 000 for the WSJT pipeline).
45    pub input_rate: u32,
46
47    /// Zero-padded forward-FFT length.
48    pub fft1_size: usize,
49
50    /// Inverse-FFT length = output sample count per call.
51    pub fft2_size: usize,
52
53    /// Tone spacing of the modulation in Hz (passed as the bandwidth unit so
54    /// the extracted band scales with the protocol).
55    pub tone_spacing_hz: f32,
56
57    /// Bins of headroom below `f0` in tone-spacing units (FT8 uses 1.5).
58    pub leading_pad_tones: f32,
59
60    /// Bins of headroom above `f0 + (ntones-1)·tone_spacing` in tone-spacing
61    /// units (FT8 uses 1.5 past tone 7 → 8.5 total).
62    pub trailing_pad_tones: f32,
63
64    /// Number of data tones; used together with `trailing_pad_tones` to
65    /// determine the upper edge of the extracted band.
66    pub ntones: u32,
67
68    /// Length of the raised-cosine taper applied to each edge (typically 101).
69    pub edge_taper_bins: usize,
70}
71
72impl DownsampleCfg {
73    #[inline]
74    fn bin_hz(&self) -> f32 {
75        self.input_rate as f32 / self.fft1_size as f32
76    }
77}
78
79/// Compute only the large forward-FFT cache. Subsequent calls to
80/// [`downsample_cached`] reuse it across all candidate frequencies, which is
81/// the expensive operation amortised by `ft8-core::decode::process_candidate`.
82#[inline]
83pub fn build_fft_cache(audio: &[i16], cfg: &DownsampleCfg) -> Vec<Complex<f32>> {
84    let mut x: Vec<Complex<f32>> = audio
85        .iter()
86        .map(|&s| Complex::new(s as f32, 0.0))
87        .chain(iter::repeat(Complex::new(0.0, 0.0)))
88        .take(cfg.fft1_size)
89        .collect();
90    let mut planner = default_planner();
91    planner.plan_forward(cfg.fft1_size).process(&mut x);
92    x
93}
94
95/// Downconvert `audio` to a complex baseband centred on `f0`.
96///
97/// Returns the decimated signal plus the forward-FFT cache so the caller can
98/// feed it to subsequent frequency-shifted calls without recomputing the
99/// 192 k-point transform.
100#[inline]
101pub fn downsample(
102    audio: &[i16],
103    f0: f32,
104    cfg: &DownsampleCfg,
105) -> (Vec<Complex<f32>>, Vec<Complex<f32>>) {
106    let cache = build_fft_cache(audio, cfg);
107    let out = downsample_cached(&cache, f0, cfg);
108    (out, cache)
109}
110
111/// Downconvert using a pre-computed forward-FFT cache.
112#[inline]
113pub fn downsample_cached(
114    fft_cache: &[Complex<f32>],
115    f0: f32,
116    cfg: &DownsampleCfg,
117) -> Vec<Complex<f32>> {
118    debug_assert_eq!(fft_cache.len(), cfg.fft1_size);
119    let mut planner = default_planner();
120
121    let df = cfg.bin_hz();
122    let baud = cfg.tone_spacing_hz;
123
124    let i0 = (f0 / df).round() as usize;
125    let ft = f0 + (cfg.ntones as f32 - 1.0 + cfg.trailing_pad_tones) * baud;
126    let fb = f0 - cfg.leading_pad_tones * baud;
127    let it = ((ft / df).round() as usize).min(cfg.fft1_size / 2);
128    let ib = ((fb / df).round() as usize).max(1);
129    let k = it.saturating_sub(ib) + 1;
130
131    let mut c1 = vec![Complex::new(0.0f32, 0.0); cfg.fft2_size];
132    for (dst, src) in c1[..k.min(cfg.fft2_size)]
133        .iter_mut()
134        .zip(fft_cache[ib..=it].iter())
135    {
136        *dst = *src;
137    }
138
139    // Raised-cosine taper on leading and trailing edges.
140    let et = cfg.edge_taper_bins;
141    if et > 1 {
142        let n = et - 1;
143        let taper: Vec<f32> = (0..et)
144            .map(|i| 0.5 * (1.0 + (i as f32 * core::f32::consts::PI / n as f32).cos()))
145            .collect();
146        for i in 0..et.min(k) {
147            c1[i] *= taper[n - i];
148        }
149        if k > et {
150            for i in 0..et {
151                c1[k - et + i] *= taper[i];
152            }
153        }
154    }
155
156    // Cyclic shift so f0 lands on DC.
157    let shift = i0.saturating_sub(ib) % cfg.fft2_size;
158    c1.rotate_left(shift);
159
160    // Inverse FFT.
161    planner.plan_inverse(cfg.fft2_size).process(&mut c1);
162
163    // Combined scale factor.
164    let fac = 1.0 / ((cfg.fft1_size as f32) * (cfg.fft2_size as f32)).sqrt();
165    for s in c1.iter_mut() {
166        *s *= fac;
167    }
168
169    c1
170}
mfsk_core/core/dsp/downsample.rs

mfsk_core/core/dsp/
downsample.rs
