flowrs-engine/crates/flowrs-core/examples/emit_frame_dsp.rs

//! Emit the frame DSP kernels' f32 output, with each sample's error bound, for the scipy
//! oracle (`tools/verify_frame_dsp.py`):
//!
//! ```sh
//! cargo run -q -p flowrs-core --example emit_frame_dsp > /tmp/frame.json
//! python3 -I tools/verify_frame_dsp.py /tmp/frame.json   # needs numpy and scipy
//! ```
//!
//! Every number is printed as the f64 equal to its f32 (`f64::from`), so scipy receives the
//! exact inputs, taps and sections the kernels ran with. Per case: the kernel, its parameters
//! and the design scipy re-derives where there is one (firwin), the input frame, the f32 output
//! (`graph::widgets::dsp_frame::kernels`, what T0, the Host program and the firmware compute) and
//! the bound `dsp_frame::bound` states: `e32`, the f32 run's rounding error against the exact
//! result, and `e64`, the f64 reference's own.

use std::collections::HashMap;
use std::fmt::Write as _;

use flowrs_core::graph::widget::Widget;
use flowrs_core::graph::widgets::dsp_frame::{
    BiquadWidget, FirWidget, Params, biquad_sections, bound, fir_taps, kernels,
};

fn signal(n: usize, seed: u64) -> Vec<f32> {
    let mut s = seed;
    (0..n)
        .map(|i| {
            s ^= s << 13;
            s ^= s >> 7;
            s ^= s << 17;
            let noise = ((s >> 11) as f64 / (1u64 << 53) as f64 * 2.0 - 1.0) as f32;
            // A tone, a slow ramp and noise: the filters' pass and stop bands both have work.
            (0.8 * (i as f32 * 0.05).sin() + 0.3 * noise + 0.5 * (i as f32 / n as f32)) * 3.3
        })
        .collect()
}

fn list(v: impl IntoIterator<Item = f64>) -> String {
    let items: Vec<String> = v.into_iter().map(|x| format!("{x:?}")).collect();
    format!("[{}]", items.join(","))
}

fn f32s(v: &[f32]) -> String {
    list(v.iter().map(|x| f64::from(*x)))
}

fn p(pairs: &[(&str, &str)]) -> HashMap<String, String> {
    pairs
        .iter()
        .map(|(k, v)| ((*k).to_string(), (*v).to_string()))
        .collect()
}

fn case(out: &mut Vec<String>, head: String, x: &[f32], y: &[f32], r: &bound::Reference) {
    let mut s = String::new();
    let _ = write!(
        s,
        "{{{head},\"x\":{},\"y\":{},\"e32\":{},\"e64\":{}}}",
        f32s(x),
        f32s(y),
        list(r.e32.iter().copied()),
        list(r.e64.iter().copied())
    );
    out.push(s);
}

fn main() {
    let mut out = Vec::new();
    let fir_defs = FirWidget.parameters();
    let bq_defs = BiquadWidget.parameters();
    for seed in [1u64, 2, 3] {
        let x = signal(256, seed);
        // FIR: scipy's firwin re-derived, every window, and written taps.
        for (window, scipy_window) in [
            ("hamming", "hamming"),
            ("hann", "hann"),
            ("blackman", "blackman"),
            ("rectangular", "boxcar"),
        ] {
            let set = p(&[
                ("num_taps", "31"),
                ("cutoff", "50"),
                ("sample_rate", "1000"),
                ("window", window),
            ]);
            let taps = fir_taps(&Params::new(&set, &fir_defs)).unwrap();
            let y = kernels::fir(&x, &taps);
            case(
                &mut out,
                format!(
                    "\"kind\":\"fir\",\"design\":{{\"numtaps\":31,\"cutoff\":50.0,\"fs\":1000.0,\"window\":\"{scipy_window}\"}},\"taps\":{}",
                    f32s(&taps)
                ),
                &x,
                &y,
                &bound::fir(&x, &taps),
            );
        }
        let written = [0.7f32, -1.3, 2.9, 0.01, -0.25];
        case(
            &mut out,
            format!("\"kind\":\"fir\",\"taps\":{}", f32s(&written)),
            &x,
            &kernels::fir(&x, &written),
            &bound::fir(&x, &written),
        );
        // Biquad cascades: every family the panel offers, both responses, several orders.
        for (family, response, order, cutoff) in [
            ("butterworth", "lowpass", "2", "50"),
            ("butterworth", "lowpass", "6", "50"),
            ("chebyshev1", "highpass", "3", "100"),
            ("chebyshev2", "lowpass", "4", "120"),
            ("elliptic", "lowpass", "6", "120"),
            ("bessel", "lowpass", "8", "20"),
        ] {
            let set = p(&[
                ("family", family),
                ("response", response),
                ("order", order),
                ("cutoff", cutoff),
            ]);
            let sos = biquad_sections(&Params::new(&set, &bq_defs)).unwrap();
            let rows: Vec<String> = sos.iter().map(|s| f32s(s)).collect();
            case(
                &mut out,
                format!(
                    "\"kind\":\"sos\",\"family\":\"{family} {response} {order}\",\"sos\":[{}]",
                    rows.join(",")
                ),
                &x,
                &kernels::biquads(&x, &sos),
                &bound::biquads(&x, &sos),
            );
        }
        // Moving statistics.
        for stat in ["mean", "rms", "min", "max", "variance"] {
            for w in [1usize, 7, 64] {
                case(
                    &mut out,
                    format!("\"kind\":\"moving\",\"stat\":\"{stat}\",\"window\":{w}"),
                    &x,
                    &kernels::moving(&x, w, stat),
                    &bound::moving(&x, w, stat),
                );
            }
        }
        // Decimation is exact.
        for f in [2usize, 3, 16] {
            let y = kernels::decimate(&x, f);
            let zero = bound::Reference {
                exact: y.iter().map(|v| f64::from(*v)).collect(),
                e32: vec![0.0; y.len()],
                e64: vec![0.0; y.len()],
            };
            case(
                &mut out,
                format!("\"kind\":\"decimate\",\"factor\":{f}"),
                &x,
                &y,
                &zero,
            );
        }
    }
    println!("{{\"cases\":[{}]}}", out.join(",\n"));
}