//! 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 { 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) -> String { let items: Vec = 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 { pairs .iter() .map(|(k, v)| ((*k).to_string(), (*v).to_string())) .collect() } fn case(out: &mut Vec, 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 = 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")); }