Frequency spectrum

use fft.qweruse fft_spectrum, fft_peaks from fft      // Quillo 0.6: only these

rustfft = "6"

// A guitar tuner's back end and a vibration check: samples in, frequencies out.// (Here the signals are made up: an A string at 110 Hz with its overtone, and a fan running at 25 Hz.)////   quillo run fft.qwe      then open http://localhost:8080/tuneruse "fft.qwer"fn tone(freqs, rate, n) = range(n) | map((i) => sum(freqs | map((f) => f[1] * math.sin(2.0 * math.pi * f[0] * i / rate))))route GET /tuner {    let samples = tone([[110.0, 0.8], [220.0, 0.3]], 8000.0, 4096)    let p = fft_peaks(samples, 8000.0, 2)?    { note_hz: round(p[0].frequency, 0), overtone_hz: round(p[1].frequency, 0), level: fft_level(samples)? }}route GET /fan {    let samples = tone([[25.0, 1.0]], 1000.0, 1000)    fft_spectrum(samples, 1000.0, false, 3)?}

//! The frequency spectrum of a signal: samples from a microphone, an accelerometer, a sensor or a price series turned into frequencies and their strength, the dominant frequency found, and a signal's loudness. For tuners, vibration monitoring and seasonality.//! version: 1.0.0//!//! [dependencies]//! rustfft = "6"use rustfft::num_complex::Complex;use serde_json::{json, Value};const MAX_SAMPLES: usize = 1 << 22;/// The one-sided amplitude spectrum, with a Hann window when asked: (frequency step, amplitudes for bins 0..=n/2).fn spectrum(samples: &[f64], rate: f64, window: bool) -> Result<(f64, Vec<f64>), String> {    let n = samples.len();    if n < 2 {        return Err("give at least 2 samples".into());    }    if n > MAX_SAMPLES {        return Err("at most 4 million samples".into());    }    if !(rate > 0.0 && rate.is_finite()) {        return Err("the sample rate is a number above 0 (samples per second)".into());    }    if samples.iter().any(|x| !x.is_finite()) {        return Err("a sample is not a finite number".into());    }    let mean = samples.iter().sum::<f64>() / n as f64;    let w: Vec<f64> = (0..n).map(|i| if window { 0.5 - 0.5 * (2.0 * std::f64::consts::PI * i as f64 / (n as f64 - 1.0)).cos() } else { 1.0 }).collect();    let gain: f64 = w.iter().sum::<f64>() / n as f64;    let mut buf: Vec<Complex<f64>> = samples.iter().zip(&w).map(|(x, w)| Complex::new((x - mean) * w, 0.0)).collect();    rustfft::FftPlanner::new().plan_fft_forward(n).process(&mut buf);    let half = n / 2;    let amps: Vec<f64> = (0..=half)        .map(|k| {            let a = buf[k].norm() / (n as f64 * gain);            if k == 0 || (n % 2 == 0 && k == half) { a } else { 2.0 * a }        })        .collect();    let mut amps = amps;    amps[0] = mean.abs();    Ok((rate / n as f64, amps))}fn r(x: f64) -> f64 {    (x * 1e6).round() / 1e6}/// The spectrum of `samples` taken `sample_rate` times a second: `[{ frequency, amplitude }]` from 0 Hz up to half the rate, in steps of rate / count./// The mean is taken out first (it is bin 0's amplitude); a Hann window is applied when `window` is true, which is best for real recordings. At most `max_bins` bins (2 to 100,000), the strongest kept, in frequency order.pub fn fft_spectrum(samples: Vec<f64>, sample_rate: f64, window: bool, max_bins: i64) -> Result<Vec<Value>, String> {    let (step, amps) = spectrum(&samples, sample_rate, window)?;    let mut idx: Vec<usize> = (0..amps.len()).collect();    let keep = (max_bins.clamp(2, 100_000) as usize).min(amps.len());    if keep < amps.len() {        idx.sort_by(|a, b| amps[*b].partial_cmp(&amps[*a]).unwrap_or(std::cmp::Ordering::Equal));        idx.truncate(keep);        idx.sort();    }    Ok(idx.into_iter().map(|k| json!({ "frequency": r(k as f64 * step), "amplitude": r(amps[k]) })).collect())}/// The strongest frequencies (not counting 0 Hz): `[{ frequency, amplitude }]`, the `count` highest peaks (1 to 100), strongest first, each refined between bins for a closer frequency./// For a tuner, the note being played; for a machine, its vibration; for daily sales, a weekly cycle (sample_rate 1 per day gives frequencies per day: 1/7 ≈ 0.143).pub fn fft_peaks(samples: Vec<f64>, sample_rate: f64, count: i64) -> Result<Vec<Value>, String> {    let (step, a) = spectrum(&samples, sample_rate, true)?;    let mut peaks: Vec<(f64, f64)> = Vec::new();    for k in 1..a.len() {        let left = a[k - 1];        let right = if k + 1 < a.len() { a[k + 1] } else { 0.0 };        if a[k] > left && a[k] >= right && a[k] > 1e-12 {            // Quadratic interpolation on the log amplitudes puts the peak between bins.            let (l, c, rr) = (left.max(1e-300).ln(), a[k].ln(), right.max(1e-300).ln());            let d = if k + 1 < a.len() { 0.5 * (l - rr) / (l - 2.0 * c + rr) } else { 0.0 };            let d = if d.is_finite() { d.clamp(-0.5, 0.5) } else { 0.0 };            peaks.push(((k as f64 + d) * step, a[k]));        }    }    peaks.sort_by(|x, y| y.1.partial_cmp(&x.1).unwrap_or(std::cmp::Ordering::Equal));    Ok(peaks.into_iter().take(count.clamp(1, 100) as usize).map(|(f, a)| json!({ "frequency": (f * 1000.0).round() / 1000.0, "amplitude": r(a) })).collect())}/// How strong a signal is: `{ rms, peak, mean, db }`, `db` being the RMS in decibels against full scale 1.0 (for audio samples from -1 to 1).pub fn fft_level(samples: Vec<f64>) -> Result<Value, String> {    if samples.is_empty() || samples.len() > MAX_SAMPLES {        return Err("give 1 to 4 million samples".into());    }    let n = samples.len() as f64;    let rms = (samples.iter().map(|x| x * x).sum::<f64>() / n).sqrt();    let peak = samples.iter().fold(0.0f64, |m, x| m.max(x.abs()));    Ok(json!({ "rms": r(rms), "peak": r(peak), "mean": r(samples.iter().sum::<f64>() / n), "db": if rms > 0.0 { json!((20.0 * rms.log10() * 100.0).round() / 100.0) } else { Value::Null } }))}