feat(dsp): spectrogram + AEC measurement toolkit (specview)
Add a dependency-free signal-analysis toolkit and a `specview` dev CLI to evaluate audio (especially echo cancellation) with objective numbers and a terminal spectrogram instead of ear alone. - src/dsp/: hand-written radix-2 FFT, Hann window, STFT, seeded test-signal generators (sine/log-sweep/white/pink/impulse), metrics (RMS/dBFS/peak/ERLE/ per-band energy), minimal WAV read+write, and a 24-bit-ANSI half-block spectrogram renderer (magma colormap, freq/time axes, dB legend, ASCII fallback). Pure layers have no I/O; only `wav` touches the filesystem. - src/bin/specview.rs: `gen` (conjure a test signal -> WAV), `show` (spectrogram + per-band energy summary), `erle` (broadband + per-band echo-return-loss between a before/after pair). - 27 unit tests (FFT correctness, ERLE landmarks, WAV round-trip, render shape). Verified end-to-end: log sweep renders as the expected exponential curve; a 20 dB-quieter copy reads +20.0 dB ERLE broadband and per band. Measurement substrate for upcoming AEC refinement (no shipped-path changes). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 4.8
@@ -15,6 +15,10 @@ path = "src/main.rs"
|
||||
name = "test_net"
|
||||
path = "src/bin/test_net.rs"
|
||||
|
||||
[[bin]]
|
||||
name = "specview"
|
||||
path = "src/bin/specview.rs"
|
||||
|
||||
[dependencies]
|
||||
anyhow = "1.0.102"
|
||||
async-trait = "0.1.89"
|
||||
|
||||
@@ -0,0 +1,258 @@
|
||||
//! `specview` — a developer CLI for the `dsp` toolkit: conjure test signals,
|
||||
//! draw spectrograms in the terminal, and measure echo-cancellation quality
|
||||
//! (ERLE) on WAV files. Not shipped in the GUI; a measurement aid.
|
||||
//!
|
||||
//! Usage:
|
||||
//! specview show <file.wav> [--width N] [--height N] [--max-hz HZ] [--floor DB] [--ascii]
|
||||
//! specview gen <kind> -o out.wav [--secs S] [--amp A] [--freq F] [--f0 F] [--f1 F] [--seed N]
|
||||
//! kinds: sine | sweep | white | pink | impulse | silence
|
||||
//! specview erle <before.wav> <after.wav> [--show]
|
||||
//!
|
||||
//! Examples:
|
||||
//! specview gen sweep -o sweep.wav --secs 3 --f0 100 --f1 12000
|
||||
//! specview show ~/peerspeak-recordings/2026-06-14/peerspeak-*.wav
|
||||
//! specview erle echo_only.wav after_aec.wav --show
|
||||
|
||||
use std::collections::HashMap;
|
||||
use std::path::Path;
|
||||
use std::process::ExitCode;
|
||||
|
||||
use peerspeak::dsp::{generators, metrics, render, stft, wav};
|
||||
|
||||
const SAMPLE_RATE: u32 = 48_000;
|
||||
|
||||
fn main() -> ExitCode {
|
||||
let args: Vec<String> = std::env::args().skip(1).collect();
|
||||
let Some(cmd) = args.first() else {
|
||||
eprintln!("{}", USAGE);
|
||||
return ExitCode::FAILURE;
|
||||
};
|
||||
|
||||
let rest = &args[1..];
|
||||
let result = match cmd.as_str() {
|
||||
"show" => cmd_show(rest),
|
||||
"gen" => cmd_gen(rest),
|
||||
"erle" => cmd_erle(rest),
|
||||
"-h" | "--help" | "help" => {
|
||||
println!("{}", USAGE);
|
||||
Ok(())
|
||||
}
|
||||
other => Err(format!("unknown command {other:?}\n\n{USAGE}")),
|
||||
};
|
||||
|
||||
match result {
|
||||
Ok(()) => ExitCode::SUCCESS,
|
||||
Err(e) => {
|
||||
eprintln!("error: {e}");
|
||||
ExitCode::FAILURE
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const USAGE: &str = "specview — spectrogram + AEC measurement tool\n\n\
|
||||
specview show <file.wav> [--width N] [--height N] [--max-hz HZ] [--floor DB] [--ascii]\n\
|
||||
specview gen <kind> -o out.wav [--secs S] [--amp A] [--freq F] [--f0 F] [--f1 F] [--seed N]\n\
|
||||
kinds: sine sweep white pink impulse silence\n\
|
||||
specview erle <before.wav> <after.wav> [--show]";
|
||||
|
||||
/// Renders a WAV as a spectrogram plus a numeric summary.
|
||||
fn cmd_show(args: &[String]) -> Result<(), String> {
|
||||
let (positional, flags) = parse_args(args);
|
||||
let path = positional.first().ok_or("show needs a <file.wav>")?;
|
||||
let data = wav::read(Path::new(path))?;
|
||||
|
||||
let opts = render::RenderOpts {
|
||||
width: flags.usize_or("width", 100),
|
||||
height: flags.usize_or("height", 30),
|
||||
max_hz: flags.f32_or("max-hz", 12_000.0),
|
||||
floor_db: flags.f32_or("floor", -80.0),
|
||||
ascii: flags.present("ascii"),
|
||||
};
|
||||
let spec = stft::analyze(&data.samples, data.sample_rate, 2048, 512);
|
||||
print!("{}", render::render(&spec, &opts));
|
||||
print_summary(path, &data.samples, &spec);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Generates a test signal and writes it to a WAV file.
|
||||
fn cmd_gen(args: &[String]) -> Result<(), String> {
|
||||
let (positional, flags) = parse_args(args);
|
||||
let kind = positional.first().ok_or("gen needs a <kind>")?;
|
||||
let out = flags
|
||||
.get("o")
|
||||
.or_else(|| flags.get("out"))
|
||||
.ok_or("gen needs -o <out.wav>")?;
|
||||
|
||||
let secs = flags.f32_or("secs", 2.0);
|
||||
let amp = flags.f32_or("amp", 0.8);
|
||||
let len = (secs * SAMPLE_RATE as f32) as usize;
|
||||
let seed = flags.u64_or("seed", 1);
|
||||
|
||||
let samples = match kind.as_str() {
|
||||
"sine" => generators::sine(flags.f32_or("freq", 1000.0), amp, len, SAMPLE_RATE),
|
||||
"sweep" => generators::log_sweep(
|
||||
flags.f32_or("f0", 100.0),
|
||||
flags.f32_or("f1", 12_000.0),
|
||||
amp,
|
||||
len,
|
||||
SAMPLE_RATE,
|
||||
),
|
||||
"white" => generators::white_noise(amp, len, seed),
|
||||
"pink" => generators::pink_noise(amp, len, seed),
|
||||
"impulse" => generators::impulse(amp, len),
|
||||
"silence" => generators::silence(len),
|
||||
other => return Err(format!("unknown kind {other:?} (sine sweep white pink impulse silence)")),
|
||||
};
|
||||
|
||||
wav::write(Path::new(out), &samples, SAMPLE_RATE)?;
|
||||
println!(
|
||||
"wrote {} — {} ({:.2}s, {} samples, rms {:.1} dBFS)",
|
||||
out,
|
||||
kind,
|
||||
secs,
|
||||
samples.len(),
|
||||
metrics::dbfs(metrics::rms(&samples)),
|
||||
);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Measures ERLE between a "before" (echo present) and "after" (post-AEC) file —
|
||||
/// broadband and per voice band — so you can see how much echo was removed and
|
||||
/// in which frequency range any residual lives.
|
||||
fn cmd_erle(args: &[String]) -> Result<(), String> {
|
||||
let (positional, flags) = parse_args(args);
|
||||
let before = positional.first().ok_or("erle needs <before.wav> <after.wav>")?;
|
||||
let after = positional.get(1).ok_or("erle needs <before.wav> <after.wav>")?;
|
||||
|
||||
let b = wav::read(Path::new(before))?;
|
||||
let a = wav::read(Path::new(after))?;
|
||||
// Compare over the overlapping length.
|
||||
let n = b.samples.len().min(a.samples.len());
|
||||
let bs = &b.samples[..n];
|
||||
let as_ = &a.samples[..n];
|
||||
|
||||
let broadband = metrics::erle(bs, as_);
|
||||
println!("ERLE (broadband): {broadband:+.1} dB (higher = more echo removed)");
|
||||
println!(
|
||||
" before: {:.1} dBFS rms after: {:.1} dBFS rms",
|
||||
metrics::dbfs(metrics::rms(bs)),
|
||||
metrics::dbfs(metrics::rms(as_)),
|
||||
);
|
||||
|
||||
// Per-band ERLE via averaged STFT magnitude energy.
|
||||
let sb = stft::analyze(bs, b.sample_rate, 2048, 512);
|
||||
let sa = stft::analyze(as_, a.sample_rate, 2048, 512);
|
||||
let avg_b = average_frame(&sb);
|
||||
let avg_a = average_frame(&sa);
|
||||
println!(" per band:");
|
||||
for band in metrics::VOICE_BANDS {
|
||||
let eb = metrics::band_energy(&avg_b, |bin| sb.bin_hz(bin), band.low_hz, band.high_hz);
|
||||
let ea = metrics::band_energy(&avg_a, |bin| sa.bin_hz(bin), band.low_hz, band.high_hz);
|
||||
let db = if eb > 1e-12 && ea > 1e-12 {
|
||||
10.0 * (eb / ea).log10()
|
||||
} else if eb > 1e-12 {
|
||||
120.0
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
println!(" {:<18} {:+6.1} dB", band.label, db);
|
||||
}
|
||||
|
||||
if flags.present("show") {
|
||||
println!("\n--- before ---");
|
||||
print!("{}", render::render(&sb, &render::RenderOpts::default()));
|
||||
println!("\n--- after ---");
|
||||
print!("{}", render::render(&sa, &render::RenderOpts::default()));
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Prints RMS/peak/dBFS plus the per-band energy distribution of a signal.
|
||||
fn print_summary(path: &str, samples: &[f32], spec: &stft::Spectrogram) {
|
||||
println!(
|
||||
"\n{}\n duration {:.2}s rms {:.1} dBFS peak {:.1} dBFS",
|
||||
path,
|
||||
samples.len() as f32 / spec.sample_rate as f32,
|
||||
metrics::dbfs(metrics::rms(samples)),
|
||||
metrics::dbfs(metrics::peak(samples)),
|
||||
);
|
||||
let avg = average_frame(spec);
|
||||
let total: f32 = avg.iter().map(|m| m * m).sum::<f32>().max(1e-12);
|
||||
println!(" energy by band:");
|
||||
for band in metrics::VOICE_BANDS {
|
||||
let e = metrics::band_energy(&avg, |bin| spec.bin_hz(bin), band.low_hz, band.high_hz);
|
||||
let pct = 100.0 * e / total;
|
||||
let bar = "█".repeat((pct / 2.5).round() as usize);
|
||||
println!(" {:<18} {:5.1}% {}", band.label, pct, bar);
|
||||
}
|
||||
}
|
||||
|
||||
/// Mean magnitude spectrum across all STFT frames (the time-averaged spectrum).
|
||||
fn average_frame(spec: &stft::Spectrogram) -> Vec<f32> {
|
||||
let bins = spec.bins();
|
||||
let mut avg = vec![0.0f32; bins];
|
||||
if spec.frames.is_empty() {
|
||||
return avg;
|
||||
}
|
||||
for frame in &spec.frames {
|
||||
for (bin, a) in avg.iter_mut().enumerate() {
|
||||
*a += frame[bin];
|
||||
}
|
||||
}
|
||||
let n = spec.frames.len() as f32;
|
||||
for a in &mut avg {
|
||||
*a /= n;
|
||||
}
|
||||
avg
|
||||
}
|
||||
|
||||
/// A trivial `--flag value` / `--flag` (boolean) parser. Positional args are
|
||||
/// everything not consumed as a flag value. Keeps the binary dependency-free.
|
||||
struct Flags {
|
||||
map: HashMap<String, String>,
|
||||
bools: Vec<String>,
|
||||
}
|
||||
|
||||
impl Flags {
|
||||
fn get(&self, key: &str) -> Option<&String> {
|
||||
self.map.get(key)
|
||||
}
|
||||
fn present(&self, key: &str) -> bool {
|
||||
self.bools.iter().any(|b| b == key) || self.map.contains_key(key)
|
||||
}
|
||||
fn f32_or(&self, key: &str, default: f32) -> f32 {
|
||||
self.map.get(key).and_then(|v| v.parse().ok()).unwrap_or(default)
|
||||
}
|
||||
fn usize_or(&self, key: &str, default: usize) -> usize {
|
||||
self.map.get(key).and_then(|v| v.parse().ok()).unwrap_or(default)
|
||||
}
|
||||
fn u64_or(&self, key: &str, default: u64) -> u64 {
|
||||
self.map.get(key).and_then(|v| v.parse().ok()).unwrap_or(default)
|
||||
}
|
||||
}
|
||||
|
||||
/// Splits args into positionals and flags. A token like `--width`/`-o` followed
|
||||
/// by a non-flag token captures that token as its value; otherwise it's a bool.
|
||||
fn parse_args(args: &[String]) -> (Vec<String>, Flags) {
|
||||
let mut positional = Vec::new();
|
||||
let mut map = HashMap::new();
|
||||
let mut bools = Vec::new();
|
||||
let mut i = 0;
|
||||
while i < args.len() {
|
||||
let tok = &args[i];
|
||||
if let Some(key) = tok.strip_prefix("--").or_else(|| tok.strip_prefix('-')) {
|
||||
let next_is_value = args.get(i + 1).is_some_and(|n| !n.starts_with('-'));
|
||||
if next_is_value {
|
||||
map.insert(key.to_string(), args[i + 1].clone());
|
||||
i += 2;
|
||||
} else {
|
||||
bools.push(key.to_string());
|
||||
i += 1;
|
||||
}
|
||||
} else {
|
||||
positional.push(tok.clone());
|
||||
i += 1;
|
||||
}
|
||||
}
|
||||
(positional, Flags { map, bools })
|
||||
}
|
||||
+186
@@ -0,0 +1,186 @@
|
||||
//! A small, dependency-free radix-2 Cooley-Tukey FFT.
|
||||
//!
|
||||
//! We hand-roll this rather than pull in `rustfft`/`realfft` because (a) the
|
||||
//! whole DSP toolkit is a developer/measurement aid, not a hot real-time path,
|
||||
//! and (b) it keeps the supply-chain surface at zero new crates. The transform
|
||||
//! is the textbook iterative in-place algorithm (bit-reversal permutation +
|
||||
//! log2(N) butterfly stages), computed in `f64` for headroom even though the
|
||||
//! audio it analyses is `f32`.
|
||||
//!
|
||||
//! Only power-of-two lengths are supported; the STFT layer always pads frames
|
||||
//! up to a power of two before calling in.
|
||||
|
||||
use std::f64::consts::PI;
|
||||
|
||||
/// A minimal complex number for the transform. Kept local (rather than pulling a
|
||||
/// `num-complex` dependency) since the FFT is the only thing that needs it.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct Complex {
|
||||
pub re: f64,
|
||||
pub im: f64,
|
||||
}
|
||||
|
||||
impl Complex {
|
||||
pub const ZERO: Complex = Complex { re: 0.0, im: 0.0 };
|
||||
|
||||
pub fn new(re: f64, im: f64) -> Self {
|
||||
Complex { re, im }
|
||||
}
|
||||
|
||||
/// Magnitude `sqrt(re^2 + im^2)`.
|
||||
pub fn magnitude(self) -> f64 {
|
||||
self.re.hypot(self.im)
|
||||
}
|
||||
|
||||
fn add(self, o: Complex) -> Complex {
|
||||
Complex::new(self.re + o.re, self.im + o.im)
|
||||
}
|
||||
|
||||
fn sub(self, o: Complex) -> Complex {
|
||||
Complex::new(self.re - o.re, self.im - o.im)
|
||||
}
|
||||
|
||||
fn mul(self, o: Complex) -> Complex {
|
||||
Complex::new(
|
||||
self.re * o.re - self.im * o.im,
|
||||
self.re * o.im + self.im * o.re,
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
/// In-place forward FFT. `buf.len()` must be a power of two.
|
||||
///
|
||||
/// Uses the standard sign convention `X[k] = sum_n x[n] * exp(-2πi·kn/N)`.
|
||||
pub fn fft(buf: &mut [Complex]) {
|
||||
transform(buf, false);
|
||||
}
|
||||
|
||||
/// In-place inverse FFT (normalized by `1/N`), the exact inverse of [`fft`].
|
||||
pub fn ifft(buf: &mut [Complex]) {
|
||||
transform(buf, true);
|
||||
let n = buf.len() as f64;
|
||||
for c in buf.iter_mut() {
|
||||
c.re /= n;
|
||||
c.im /= n;
|
||||
}
|
||||
}
|
||||
|
||||
fn transform(buf: &mut [Complex], inverse: bool) {
|
||||
let n = buf.len();
|
||||
assert!(n.is_power_of_two(), "FFT length {n} must be a power of two");
|
||||
if n <= 1 {
|
||||
return;
|
||||
}
|
||||
|
||||
// Bit-reversal permutation: reorder so the iterative butterflies can run
|
||||
// bottom-up in place.
|
||||
let mut j = 0usize;
|
||||
for i in 1..n {
|
||||
let mut bit = n >> 1;
|
||||
while j & bit != 0 {
|
||||
j ^= bit;
|
||||
bit >>= 1;
|
||||
}
|
||||
j ^= bit;
|
||||
if i < j {
|
||||
buf.swap(i, j);
|
||||
}
|
||||
}
|
||||
|
||||
// Butterfly stages: combine length-`len` DFTs from length-`len/2` halves,
|
||||
// doubling `len` each pass.
|
||||
let sign = if inverse { 1.0 } else { -1.0 };
|
||||
let mut len = 2;
|
||||
while len <= n {
|
||||
let ang = sign * 2.0 * PI / len as f64;
|
||||
let wlen = Complex::new(ang.cos(), ang.sin());
|
||||
let mut i = 0;
|
||||
while i < n {
|
||||
let mut w = Complex::new(1.0, 0.0);
|
||||
for k in 0..len / 2 {
|
||||
let u = buf[i + k];
|
||||
let v = buf[i + k + len / 2].mul(w);
|
||||
buf[i + k] = u.add(v);
|
||||
buf[i + k + len / 2] = u.sub(v);
|
||||
w = w.mul(wlen);
|
||||
}
|
||||
i += len;
|
||||
}
|
||||
len <<= 1;
|
||||
}
|
||||
}
|
||||
|
||||
/// Forward FFT of a real signal, returning the **one-sided** magnitude spectrum:
|
||||
/// bins `0..=N/2` (DC through Nyquist), where `N` is the next power of two ≥
|
||||
/// `samples.len()`. The input is zero-padded up to `N`.
|
||||
///
|
||||
/// Magnitudes are raw (un-normalized) linear amplitudes; callers convert to dB
|
||||
/// or normalize as needed.
|
||||
pub fn real_magnitude_spectrum(samples: &[f32]) -> Vec<f32> {
|
||||
let n = samples.len().next_power_of_two().max(2);
|
||||
let mut buf = vec![Complex::ZERO; n];
|
||||
for (i, &s) in samples.iter().enumerate() {
|
||||
buf[i].re = s as f64;
|
||||
}
|
||||
fft(&mut buf);
|
||||
buf[..=n / 2].iter().map(|c| c.magnitude() as f32).collect()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn approx(a: f64, b: f64, eps: f64) -> bool {
|
||||
(a - b).abs() <= eps
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn impulse_transforms_to_flat_spectrum() {
|
||||
// FFT of a unit impulse at n=0 is all-ones (flat spectrum).
|
||||
let mut buf = vec![Complex::ZERO; 8];
|
||||
buf[0] = Complex::new(1.0, 0.0);
|
||||
fft(&mut buf);
|
||||
for c in &buf {
|
||||
assert!(approx(c.magnitude(), 1.0, 1e-9), "expected flat 1.0, got {c:?}");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn single_bin_sine_peaks_in_that_bin() {
|
||||
// A cosine at exactly bin k=2 over N=16 should put all energy in bin 2
|
||||
// (and its mirror N-2). Check the one-sided spectrum peaks at bin 2.
|
||||
let n = 16;
|
||||
let k = 2;
|
||||
let samples: Vec<f32> = (0..n)
|
||||
.map(|i| (2.0 * PI * k as f64 * i as f64 / n as f64).cos() as f32)
|
||||
.collect();
|
||||
let mag = real_magnitude_spectrum(&samples);
|
||||
let peak_bin = mag
|
||||
.iter()
|
||||
.enumerate()
|
||||
.max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
|
||||
.unwrap()
|
||||
.0;
|
||||
assert_eq!(peak_bin, k, "energy should land in bin {k}, got {peak_bin}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn ifft_inverts_fft() {
|
||||
let original: Vec<Complex> = (0..32)
|
||||
.map(|i| Complex::new((i as f64 * 0.3).sin(), (i as f64 * 0.1).cos()))
|
||||
.collect();
|
||||
let mut buf = original.clone();
|
||||
fft(&mut buf);
|
||||
ifft(&mut buf);
|
||||
for (a, b) in original.iter().zip(&buf) {
|
||||
assert!(approx(a.re, b.re, 1e-9) && approx(a.im, b.im, 1e-9));
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[should_panic(expected = "power of two")]
|
||||
fn non_power_of_two_panics() {
|
||||
let mut buf = vec![Complex::ZERO; 6];
|
||||
fft(&mut buf);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,156 @@
|
||||
//! Reproducible test-signal generators — the "conjured noises" the measurement
|
||||
//! harness drives AEC/analysis with. Everything is deterministic: the noise
|
||||
//! sources take an explicit seed so a run is byte-for-byte repeatable, which is
|
||||
//! what makes A/B comparisons (old build vs new, tuning X vs Y) meaningful.
|
||||
//!
|
||||
//! All generators emit `f32` samples nominally in `[-amp, amp]` at the given
|
||||
//! sample rate. Convert to `i16` at the I/O edge.
|
||||
|
||||
use std::f64::consts::PI;
|
||||
|
||||
/// A tiny deterministic PRNG (xorshift64*) so the noise generators need no
|
||||
/// external crate and produce identical output for a given seed across runs.
|
||||
struct XorShift64 {
|
||||
state: u64,
|
||||
}
|
||||
|
||||
impl XorShift64 {
|
||||
fn new(seed: u64) -> Self {
|
||||
// Avoid the all-zero fixed point.
|
||||
XorShift64 {
|
||||
state: seed.wrapping_mul(0x9E37_79B9_7F4A_7C15) | 1,
|
||||
}
|
||||
}
|
||||
|
||||
fn next_u64(&mut self) -> u64 {
|
||||
let mut x = self.state;
|
||||
x ^= x >> 12;
|
||||
x ^= x << 25;
|
||||
x ^= x >> 27;
|
||||
self.state = x;
|
||||
x.wrapping_mul(0x2545_F491_4F6C_DD1D)
|
||||
}
|
||||
|
||||
/// Uniform `f32` in `[-1.0, 1.0)`.
|
||||
fn next_bipolar(&mut self) -> f32 {
|
||||
// Top 24 bits -> [0,1), then map to [-1,1).
|
||||
let u = (self.next_u64() >> 40) as f32 / (1u32 << 24) as f32;
|
||||
u * 2.0 - 1.0
|
||||
}
|
||||
}
|
||||
|
||||
/// `len` samples of silence.
|
||||
pub fn silence(len: usize) -> Vec<f32> {
|
||||
vec![0.0; len]
|
||||
}
|
||||
|
||||
/// A pure sine tone at `freq` Hz.
|
||||
pub fn sine(freq: f32, amp: f32, len: usize, sample_rate: u32) -> Vec<f32> {
|
||||
(0..len)
|
||||
.map(|i| amp * (2.0 * PI * freq as f64 * i as f64 / sample_rate as f64).sin() as f32)
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// A unit impulse: `amp` at sample 0, silence after. Its flat spectrum makes it
|
||||
/// the natural probe for an echo path's impulse response.
|
||||
pub fn impulse(amp: f32, len: usize) -> Vec<f32> {
|
||||
let mut v = vec![0.0; len];
|
||||
if len > 0 {
|
||||
v[0] = amp;
|
||||
}
|
||||
v
|
||||
}
|
||||
|
||||
/// A logarithmic (exponential) sine sweep from `f0` to `f1` Hz over the buffer —
|
||||
/// the standard excitation for measuring a system's frequency response, since it
|
||||
/// spends equal time per octave.
|
||||
pub fn log_sweep(f0: f32, f1: f32, amp: f32, len: usize, sample_rate: u32) -> Vec<f32> {
|
||||
if len == 0 {
|
||||
return Vec::new();
|
||||
}
|
||||
let t_total = len as f64 / sample_rate as f64;
|
||||
let (f0, f1) = (f0 as f64, f1 as f64);
|
||||
let k = (f1 / f0).ln();
|
||||
(0..len)
|
||||
.map(|i| {
|
||||
let t = i as f64 / sample_rate as f64;
|
||||
// Instantaneous-phase integral of an exponential chirp.
|
||||
let phase = 2.0 * PI * f0 * t_total / k * ((k * t / t_total).exp() - 1.0);
|
||||
amp * phase.sin() as f32
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Uniform white noise in `[-amp, amp]`, seeded for reproducibility.
|
||||
pub fn white_noise(amp: f32, len: usize, seed: u64) -> Vec<f32> {
|
||||
let mut rng = XorShift64::new(seed);
|
||||
(0..len).map(|_| amp * rng.next_bipolar()).collect()
|
||||
}
|
||||
|
||||
/// Pink (1/f) noise via the Voss-McCartney algorithm — perceptually flatter than
|
||||
/// white noise and a closer stand-in for room/voice spectra. Seeded.
|
||||
pub fn pink_noise(amp: f32, len: usize, seed: u64) -> Vec<f32> {
|
||||
let mut rng = XorShift64::new(seed);
|
||||
const ROWS: usize = 16;
|
||||
let mut rows = [0.0f32; ROWS];
|
||||
let mut running = 0.0f32;
|
||||
let mut out = Vec::with_capacity(len);
|
||||
for i in 0..len {
|
||||
// Each row updates half as often as the previous: the row to refresh is
|
||||
// the index of the lowest set bit of the sample counter.
|
||||
let n = i + 1;
|
||||
let row = (n & n.wrapping_neg()).trailing_zeros() as usize;
|
||||
if row < ROWS {
|
||||
running -= rows[row];
|
||||
rows[row] = rng.next_bipolar();
|
||||
running += rows[row];
|
||||
}
|
||||
// Normalize: sum of ROWS unit-bipolar rows spans roughly [-ROWS, ROWS].
|
||||
out.push(amp * running / ROWS as f32);
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::dsp::metrics::rms;
|
||||
|
||||
#[test]
|
||||
fn sine_has_expected_length_and_range() {
|
||||
let s = sine(1000.0, 0.5, 480, 48_000);
|
||||
assert_eq!(s.len(), 480);
|
||||
assert!(s.iter().all(|&x| x.abs() <= 0.5001));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn impulse_is_single_spike() {
|
||||
let s = impulse(1.0, 10);
|
||||
assert_eq!(s[0], 1.0);
|
||||
assert!(s[1..].iter().all(|&x| x == 0.0));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn white_noise_is_deterministic_for_a_seed() {
|
||||
let a = white_noise(1.0, 256, 42);
|
||||
let b = white_noise(1.0, 256, 42);
|
||||
let c = white_noise(1.0, 256, 43);
|
||||
assert_eq!(a, b, "same seed -> identical output");
|
||||
assert_ne!(a, c, "different seed -> different output");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noise_energy_is_nonzero_and_bounded() {
|
||||
let w = white_noise(0.5, 4096, 7);
|
||||
let r = rms(&w);
|
||||
assert!(r > 0.0 && r < 0.5, "rms {r} should be in (0, amp)");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn log_sweep_spans_the_buffer() {
|
||||
let s = log_sweep(100.0, 8000.0, 1.0, 4800, 48_000);
|
||||
assert_eq!(s.len(), 4800);
|
||||
// Non-trivial energy present.
|
||||
assert!(rms(&s) > 0.1);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,132 @@
|
||||
//! Objective signal metrics — the numbers that turn "sounds better" into a
|
||||
//! measurement. The headline one for echo cancellation is **ERLE** (Echo Return
|
||||
//! Loss Enhancement): how much echo energy the canceller removed, in dB.
|
||||
|
||||
/// Root-mean-square level of a signal (linear amplitude). `0.0` for empty input.
|
||||
pub fn rms(samples: &[f32]) -> f32 {
|
||||
if samples.is_empty() {
|
||||
return 0.0;
|
||||
}
|
||||
let sum_sq: f64 = samples.iter().map(|&s| s as f64 * s as f64).sum();
|
||||
(sum_sq / samples.len() as f64).sqrt() as f32
|
||||
}
|
||||
|
||||
/// Peak absolute amplitude. `0.0` for empty input.
|
||||
pub fn peak(samples: &[f32]) -> f32 {
|
||||
samples.iter().fold(0.0, |m, &s| m.max(s.abs()))
|
||||
}
|
||||
|
||||
/// Converts a linear amplitude (e.g. an RMS or peak value, relative to full-scale
|
||||
/// `1.0`) to decibels below full scale. A floor of -120 dBFS is returned for
|
||||
/// silence so the result is always finite.
|
||||
pub fn dbfs(linear: f32) -> f32 {
|
||||
if linear <= 1e-6 {
|
||||
return -120.0;
|
||||
}
|
||||
20.0 * linear.log10()
|
||||
}
|
||||
|
||||
/// **Echo Return Loss Enhancement**, in dB: `10·log10(E[echo²] / E[residual²])`.
|
||||
///
|
||||
/// `echo` is the signal *before* cancellation (the echo the mic picked up);
|
||||
/// `residual` is what's *left after* the canceller ran. A larger number is
|
||||
/// better — e.g. +30 dB means the canceller removed 99.9% of the echo energy.
|
||||
/// Returns a +120 dB ceiling if the residual is effectively silent (perfect
|
||||
/// cancellation) and 0.0 if there was no echo energy to begin with.
|
||||
pub fn erle(echo: &[f32], residual: &[f32]) -> f32 {
|
||||
let e_echo = mean_square(echo);
|
||||
let e_res = mean_square(residual);
|
||||
if e_echo <= 1e-12 {
|
||||
return 0.0;
|
||||
}
|
||||
if e_res <= 1e-12 {
|
||||
return 120.0;
|
||||
}
|
||||
10.0 * (e_echo / e_res).log10() as f32
|
||||
}
|
||||
|
||||
/// Mean square (average energy per sample) of a signal. `0.0` for empty input.
|
||||
fn mean_square(samples: &[f32]) -> f64 {
|
||||
if samples.is_empty() {
|
||||
return 0.0;
|
||||
}
|
||||
samples.iter().map(|&s| s as f64 * s as f64).sum::<f64>() / samples.len() as f64
|
||||
}
|
||||
|
||||
/// A frequency band for per-band energy analysis, in Hz.
|
||||
pub struct Band {
|
||||
pub label: &'static str,
|
||||
pub low_hz: f32,
|
||||
pub high_hz: f32,
|
||||
}
|
||||
|
||||
/// Voice-relevant bands for spotting *where* residual echo or noise lives.
|
||||
pub const VOICE_BANDS: &[Band] = &[
|
||||
Band { label: "low (80-300)", low_hz: 80.0, high_hz: 300.0 },
|
||||
Band { label: "low-mid (300-1k)", low_hz: 300.0, high_hz: 1000.0 },
|
||||
Band { label: "mid (1k-3k)", low_hz: 1000.0, high_hz: 3000.0 },
|
||||
Band { label: "high-mid (3k-6k)", low_hz: 3000.0, high_hz: 6000.0 },
|
||||
Band { label: "high (6k-12k)", low_hz: 6000.0, high_hz: 12000.0 },
|
||||
];
|
||||
|
||||
/// Sums the linear magnitude energy within `[low_hz, high_hz)` across a single
|
||||
/// STFT magnitude frame. `bin_hz` maps a bin index to its centre frequency.
|
||||
pub fn band_energy(frame: &[f32], bin_hz: impl Fn(usize) -> f32, low_hz: f32, high_hz: f32) -> f32 {
|
||||
frame
|
||||
.iter()
|
||||
.enumerate()
|
||||
.filter(|&(bin, _)| {
|
||||
let hz = bin_hz(bin);
|
||||
hz >= low_hz && hz < high_hz
|
||||
})
|
||||
.map(|(_, &m)| m * m)
|
||||
.sum()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn rms_of_constant_is_that_constant() {
|
||||
assert!((rms(&[0.5; 100]) - 0.5).abs() < 1e-6);
|
||||
assert_eq!(rms(&[]), 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn peak_finds_max_magnitude() {
|
||||
assert_eq!(peak(&[0.1, -0.9, 0.3]), 0.9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn dbfs_landmarks() {
|
||||
assert!((dbfs(1.0) - 0.0).abs() < 1e-4, "full scale = 0 dBFS");
|
||||
assert!((dbfs(0.5) - -6.0206).abs() < 1e-3, "half = ~-6 dB");
|
||||
assert_eq!(dbfs(0.0), -120.0, "silence floors");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn erle_halving_energy_is_about_3db() {
|
||||
// residual amplitude = echo/sqrt(2) -> half the energy -> ~3.01 dB.
|
||||
let echo = vec![1.0f32; 1000];
|
||||
let residual = vec![std::f32::consts::FRAC_1_SQRT_2; 1000];
|
||||
let e = erle(&echo, &residual);
|
||||
assert!((e - 3.0103).abs() < 0.01, "expected ~3 dB, got {e}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn erle_perfect_cancellation_ceils() {
|
||||
assert_eq!(erle(&[1.0; 10], &[0.0; 10]), 120.0);
|
||||
assert_eq!(erle(&[0.0; 10], &[0.0; 10]), 0.0, "no echo -> 0");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn band_energy_selects_the_right_bins() {
|
||||
// 5-bin frame, 100 Hz per bin: bins at 0,100,200,300,400 Hz.
|
||||
let frame = [1.0, 2.0, 3.0, 4.0, 5.0];
|
||||
let bin_hz = |b: usize| b as f32 * 100.0;
|
||||
// [150,350) -> bins 200,300 Hz -> 3^2 + 4^2 = 25.
|
||||
let e = band_energy(&frame, bin_hz, 150.0, 350.0);
|
||||
assert!((e - 25.0).abs() < 1e-4, "got {e}");
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,18 @@
|
||||
//! Signal-analysis toolkit: FFT, STFT spectrograms, reproducible test-signal
|
||||
//! generators, objective metrics (RMS/dBFS/ERLE/per-band energy), a minimal WAV
|
||||
//! reader/writer, and a 24-bit-colour terminal spectrogram renderer.
|
||||
//!
|
||||
//! This is a **measurement/developer aid**, not part of the real-time audio
|
||||
//! path. It exists so audio work — especially echo-cancellation tuning — can be
|
||||
//! evaluated with reproducible numbers and a visible spectrogram instead of
|
||||
//! ear alone. The pure layers (`fft`, `window`, `stft`, `generators`, `metrics`,
|
||||
//! `render`) have no I/O; `wav` is the only edge that touches the filesystem.
|
||||
//! Driven from the `specview` binary (`src/bin/specview.rs`).
|
||||
|
||||
pub mod fft;
|
||||
pub mod generators;
|
||||
pub mod metrics;
|
||||
pub mod render;
|
||||
pub mod stft;
|
||||
pub mod wav;
|
||||
pub mod window;
|
||||
@@ -0,0 +1,260 @@
|
||||
//! Renders a [`Spectrogram`] to a string for the terminal.
|
||||
//!
|
||||
//! The trick that gives real resolution in a text grid: each character cell is
|
||||
//! the Unicode upper-half block `▀` with a 24-bit-colour **foreground** (top
|
||||
//! pixel) and **background** (bottom pixel), so one row of text shows two
|
||||
//! frequency bins. Magnitudes are mapped to dB against the loudest cell and
|
||||
//! coloured with a magma-style ramp (dark = quiet, bright = loud). An ASCII
|
||||
//! fallback ramp is available for non-truecolor / piped output.
|
||||
|
||||
use super::stft::Spectrogram;
|
||||
|
||||
/// Rendering options for [`render`].
|
||||
pub struct RenderOpts {
|
||||
/// Output width in character columns (time axis).
|
||||
pub width: usize,
|
||||
/// Output height in character rows; vertical resolution is `2×height` bins.
|
||||
pub height: usize,
|
||||
/// Dynamic-range floor in dB below the loudest cell (e.g. `-80.0`).
|
||||
pub floor_db: f32,
|
||||
/// Top of the frequency axis in Hz (clamped to Nyquist).
|
||||
pub max_hz: f32,
|
||||
/// Use a plain-ASCII brightness ramp instead of 24-bit colour blocks.
|
||||
pub ascii: bool,
|
||||
}
|
||||
|
||||
impl Default for RenderOpts {
|
||||
fn default() -> Self {
|
||||
RenderOpts {
|
||||
width: 100,
|
||||
height: 30,
|
||||
floor_db: -80.0,
|
||||
max_hz: 12_000.0,
|
||||
ascii: false,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const RESET: &str = "\x1b[0m";
|
||||
/// Brightness ramp (dark→bright) for the ASCII fallback.
|
||||
const ASCII_RAMP: &[u8] = b" .:-=+*#%@";
|
||||
|
||||
/// Renders `spec` to a multi-line string: a header, the frequency-labelled
|
||||
/// spectrogram body, a time axis, and a dB colour legend.
|
||||
pub fn render(spec: &Spectrogram, opts: &RenderOpts) -> String {
|
||||
if spec.frames.is_empty() {
|
||||
return "(empty signal — nothing to plot)\n".to_string();
|
||||
}
|
||||
|
||||
let nyquist = spec.sample_rate as f32 / 2.0;
|
||||
let max_hz = opts.max_hz.min(nyquist).max(spec.bin_hz(1));
|
||||
// Number of bins from DC up to max_hz (at least 2 so we can aggregate).
|
||||
let bins_used = (0..spec.bins())
|
||||
.take_while(|&b| spec.bin_hz(b) <= max_hz)
|
||||
.count()
|
||||
.max(2);
|
||||
|
||||
let n_frames = spec.frames.len();
|
||||
let width = opts.width.clamp(1, n_frames);
|
||||
let rows = opts.height.max(1) * 2; // 2 frequency pixels per character row
|
||||
|
||||
// Stage 1: aggregate frames into `width` time columns (mean magnitude/bin).
|
||||
let mut columns: Vec<Vec<f32>> = Vec::with_capacity(width);
|
||||
for c in 0..width {
|
||||
let f_lo = c * n_frames / width;
|
||||
let f_hi = ((c + 1) * n_frames / width).max(f_lo + 1).min(n_frames);
|
||||
let span = (f_hi - f_lo) as f32;
|
||||
let mut acc = vec![0.0f32; bins_used];
|
||||
for f in f_lo..f_hi {
|
||||
let frame = &spec.frames[f];
|
||||
for (bin, a) in acc.iter_mut().enumerate() {
|
||||
*a += frame[bin];
|
||||
}
|
||||
}
|
||||
for a in &mut acc {
|
||||
*a /= span;
|
||||
}
|
||||
columns.push(acc);
|
||||
}
|
||||
|
||||
// Stage 2: aggregate each column's bins into `rows` vertical pixels.
|
||||
// grid[col][pixel], pixel 0 = lowest freq.
|
||||
let mut grid: Vec<Vec<f32>> = Vec::with_capacity(width);
|
||||
let mut max_mag = 1e-9f32;
|
||||
for col in &columns {
|
||||
let mut pixels = vec![0.0f32; rows];
|
||||
for (p, px) in pixels.iter_mut().enumerate() {
|
||||
let b_lo = p * bins_used / rows;
|
||||
let b_hi = ((p + 1) * bins_used / rows).max(b_lo + 1).min(bins_used);
|
||||
let mean = col[b_lo..b_hi].iter().copied().sum::<f32>() / (b_hi - b_lo) as f32;
|
||||
*px = mean;
|
||||
max_mag = max_mag.max(mean);
|
||||
}
|
||||
grid.push(pixels);
|
||||
}
|
||||
|
||||
// Stage 3: draw. Char row r (0 = top) shows pixels (top=high freq, bottom).
|
||||
let mut out = String::new();
|
||||
out.push_str(&format!(
|
||||
"Spectrogram {:.0} Hz–{:.0} Hz · {:.2}s · {} frames, fft {}, hop {}\n",
|
||||
0.0,
|
||||
max_hz,
|
||||
spec.frame_time(n_frames),
|
||||
n_frames,
|
||||
spec.fft_size,
|
||||
spec.hop,
|
||||
));
|
||||
|
||||
let label_every = (opts.height.max(1) / 6).max(1);
|
||||
for r in 0..opts.height {
|
||||
let top_pixel = rows - 1 - 2 * r;
|
||||
let bot_pixel = rows.saturating_sub(2 + 2 * r);
|
||||
// Frequency label for this char row (use the top pixel's centre freq).
|
||||
let label = if r % label_every == 0 {
|
||||
let frac = top_pixel as f32 / (rows - 1) as f32;
|
||||
format!("{:>6.0}", frac * max_hz)
|
||||
} else {
|
||||
" ".to_string()
|
||||
};
|
||||
out.push_str(&label);
|
||||
out.push_str(" \u{2502}"); // " │"
|
||||
for col in &grid {
|
||||
let top = norm_db(col[top_pixel], max_mag, opts.floor_db);
|
||||
let bot = norm_db(col[bot_pixel], max_mag, opts.floor_db);
|
||||
push_cell(&mut out, top, bot, opts.ascii);
|
||||
}
|
||||
if !opts.ascii {
|
||||
out.push_str(RESET);
|
||||
}
|
||||
out.push('\n');
|
||||
}
|
||||
|
||||
// Time axis.
|
||||
out.push_str(" \u{2514}"); // " └"
|
||||
out.push_str(&"\u{2500}".repeat(width));
|
||||
out.push('\n');
|
||||
out.push_str(&format!(
|
||||
" 0.00s{:>width$.2}s\n",
|
||||
spec.frame_time(n_frames),
|
||||
width = width.saturating_sub(5),
|
||||
));
|
||||
|
||||
out.push_str(&legend(opts));
|
||||
out
|
||||
}
|
||||
|
||||
/// Maps a linear magnitude to `[0,1]` over `[floor_db, 0]` relative to `max_mag`.
|
||||
fn norm_db(mag: f32, max_mag: f32, floor_db: f32) -> f32 {
|
||||
if mag <= 1e-9 {
|
||||
return 0.0;
|
||||
}
|
||||
let db = 20.0 * (mag / max_mag).log10();
|
||||
((db - floor_db) / -floor_db).clamp(0.0, 1.0)
|
||||
}
|
||||
|
||||
/// Appends one character cell (two stacked pixels) to `out`.
|
||||
fn push_cell(out: &mut String, top: f32, bot: f32, ascii: bool) {
|
||||
if ascii {
|
||||
// One pixel per cell in ASCII mode (use the brighter of the pair).
|
||||
let v = top.max(bot);
|
||||
let idx = ((v * (ASCII_RAMP.len() - 1) as f32).round() as usize).min(ASCII_RAMP.len() - 1);
|
||||
out.push(ASCII_RAMP[idx] as char);
|
||||
} else {
|
||||
let (tr, tg, tb) = magma(top);
|
||||
let (br, bg, bb) = magma(bot);
|
||||
// fg = top pixel, bg = bottom pixel, glyph = upper half block.
|
||||
out.push_str(&format!(
|
||||
"\x1b[38;2;{tr};{tg};{tb}m\x1b[48;2;{br};{bg};{bb}m\u{2580}"
|
||||
));
|
||||
}
|
||||
}
|
||||
|
||||
/// A magma-style colormap: `t` in `[0,1]` → `(r,g,b)`. Piecewise-linear through
|
||||
/// control colours sampled from matplotlib's `magma`.
|
||||
fn magma(t: f32) -> (u8, u8, u8) {
|
||||
const STOPS: &[(f32, (f32, f32, f32))] = &[
|
||||
(0.0, (0.0, 0.0, 4.0)),
|
||||
(0.25, (80.0, 18.0, 123.0)),
|
||||
(0.5, (182.0, 54.0, 121.0)),
|
||||
(0.75, (252.0, 136.0, 97.0)),
|
||||
(1.0, (252.0, 253.0, 191.0)),
|
||||
];
|
||||
let t = t.clamp(0.0, 1.0);
|
||||
for win in STOPS.windows(2) {
|
||||
let (t0, c0) = win[0];
|
||||
let (t1, c1) = win[1];
|
||||
if t <= t1 {
|
||||
let f = if t1 > t0 { (t - t0) / (t1 - t0) } else { 0.0 };
|
||||
return (
|
||||
(c0.0 + (c1.0 - c0.0) * f) as u8,
|
||||
(c0.1 + (c1.1 - c0.1) * f) as u8,
|
||||
(c0.2 + (c1.2 - c0.2) * f) as u8,
|
||||
);
|
||||
}
|
||||
}
|
||||
(252, 253, 191)
|
||||
}
|
||||
|
||||
/// A horizontal dB colour legend strip.
|
||||
fn legend(opts: &RenderOpts) -> String {
|
||||
let mut s = format!(" dB: {:.0}", opts.floor_db);
|
||||
let steps = 32;
|
||||
if opts.ascii {
|
||||
for i in 0..steps {
|
||||
let v = i as f32 / (steps - 1) as f32;
|
||||
let idx = ((v * (ASCII_RAMP.len() - 1) as f32).round() as usize).min(ASCII_RAMP.len() - 1);
|
||||
s.push(ASCII_RAMP[idx] as char);
|
||||
}
|
||||
} else {
|
||||
s.push(' ');
|
||||
for i in 0..steps {
|
||||
let v = i as f32 / (steps - 1) as f32;
|
||||
let (r, g, b) = magma(v);
|
||||
s.push_str(&format!("\x1b[48;2;{r};{g};{b}m "));
|
||||
}
|
||||
s.push_str(RESET);
|
||||
}
|
||||
s.push_str(" 0\n");
|
||||
s
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::dsp::{generators, stft};
|
||||
|
||||
#[test]
|
||||
fn norm_db_maps_floor_and_peak() {
|
||||
// At max magnitude -> 0 dB -> top of range (1.0).
|
||||
assert!((norm_db(1.0, 1.0, -80.0) - 1.0).abs() < 1e-5);
|
||||
// Silence -> 0.0.
|
||||
assert_eq!(norm_db(0.0, 1.0, -80.0), 0.0);
|
||||
// -80 dB (mag 1e-4 of max) -> ~0.0 bottom of range.
|
||||
assert!(norm_db(1e-4, 1.0, -80.0) < 0.01);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn magma_endpoints() {
|
||||
assert_eq!(magma(0.0), (0, 0, 4));
|
||||
assert_eq!(magma(1.0), (252, 253, 191));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn render_produces_grid_of_expected_height() {
|
||||
let sig = generators::sine(2000.0, 0.8, 48_000, 48_000);
|
||||
let spec = stft::analyze(&sig, 48_000, 1024, 512);
|
||||
let opts = RenderOpts { width: 40, height: 10, ..Default::default() };
|
||||
let out = render(&spec, &opts);
|
||||
// Header + 10 body rows + time axis (2) + legend = non-trivial.
|
||||
let lines = out.lines().count();
|
||||
assert!(lines >= 13, "expected >=13 lines, got {lines}");
|
||||
assert!(out.contains("Spectrogram"));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn empty_spectrogram_renders_message() {
|
||||
let spec = stft::analyze(&[], 48_000, 512, 256);
|
||||
let out = render(&spec, &RenderOpts::default());
|
||||
assert!(out.contains("empty"));
|
||||
}
|
||||
}
|
||||
+105
@@ -0,0 +1,105 @@
|
||||
//! Short-Time Fourier Transform: slice a signal into overlapping windowed
|
||||
//! frames and take each frame's magnitude spectrum. The result is the
|
||||
//! time×frequency matrix a spectrogram draws.
|
||||
|
||||
use super::fft::real_magnitude_spectrum;
|
||||
use super::window::{apply, hann};
|
||||
|
||||
/// One STFT analysis: a sequence of magnitude frames plus the geometry needed to
|
||||
/// label axes (sample rate, fft size, hop).
|
||||
pub struct Spectrogram {
|
||||
/// `frames[t][bin]` = linear magnitude of frequency `bin` at time-step `t`.
|
||||
/// Each inner vec has `fft_size/2 + 1` bins (DC..Nyquist).
|
||||
pub frames: Vec<Vec<f32>>,
|
||||
pub sample_rate: u32,
|
||||
pub fft_size: usize,
|
||||
pub hop: usize,
|
||||
}
|
||||
|
||||
impl Spectrogram {
|
||||
/// Number of one-sided frequency bins per frame (`fft_size/2 + 1`).
|
||||
pub fn bins(&self) -> usize {
|
||||
self.fft_size / 2 + 1
|
||||
}
|
||||
|
||||
/// Centre frequency (Hz) of bin index `bin`.
|
||||
pub fn bin_hz(&self, bin: usize) -> f32 {
|
||||
bin as f32 * self.sample_rate as f32 / self.fft_size as f32
|
||||
}
|
||||
|
||||
/// Time (seconds) at the start of frame `t`.
|
||||
pub fn frame_time(&self, t: usize) -> f32 {
|
||||
(t * self.hop) as f32 / self.sample_rate as f32
|
||||
}
|
||||
}
|
||||
|
||||
/// Computes the STFT of `samples` with the given `fft_size` (rounded up to a
|
||||
/// power of two) and `hop` (frame advance in samples). A Hann window is applied
|
||||
/// to each frame. The final partial frame is zero-padded so trailing audio is
|
||||
/// not dropped.
|
||||
pub fn analyze(samples: &[f32], sample_rate: u32, fft_size: usize, hop: usize) -> Spectrogram {
|
||||
let fft_size = fft_size.next_power_of_two().max(2);
|
||||
let hop = hop.max(1);
|
||||
let window = hann(fft_size);
|
||||
|
||||
let mut frames = Vec::new();
|
||||
if !samples.is_empty() {
|
||||
let mut start = 0;
|
||||
while start < samples.len() {
|
||||
let end = (start + fft_size).min(samples.len());
|
||||
let mut frame = vec![0.0f32; fft_size];
|
||||
frame[..end - start].copy_from_slice(&samples[start..end]);
|
||||
let windowed = apply(&frame, &window);
|
||||
frames.push(real_magnitude_spectrum(&windowed));
|
||||
start += hop;
|
||||
}
|
||||
}
|
||||
|
||||
Spectrogram {
|
||||
frames,
|
||||
sample_rate,
|
||||
fft_size,
|
||||
hop,
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use std::f64::consts::PI;
|
||||
|
||||
#[test]
|
||||
fn frame_count_follows_hop() {
|
||||
// 1000 samples, hop 250 -> frames start at 0,250,500,750 = 4 frames.
|
||||
let sig = vec![0.0f32; 1000];
|
||||
let s = analyze(&sig, 48_000, 512, 250);
|
||||
assert_eq!(s.frames.len(), 4);
|
||||
assert_eq!(s.bins(), 512 / 2 + 1);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn tone_lands_in_expected_bin() {
|
||||
// A 3 kHz tone at 48 kHz with a 1024-pt FFT -> bin ≈ 3000/(48000/1024) = 64.
|
||||
let sr = 48_000;
|
||||
let freq = 3000.0;
|
||||
let sig: Vec<f32> = (0..4096)
|
||||
.map(|i| (2.0 * PI * freq * i as f64 / sr as f64).sin() as f32)
|
||||
.collect();
|
||||
let s = analyze(&sig, sr, 1024, 512);
|
||||
let mid = &s.frames[s.frames.len() / 2];
|
||||
let peak_bin = mid
|
||||
.iter()
|
||||
.enumerate()
|
||||
.max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
|
||||
.unwrap()
|
||||
.0;
|
||||
let peak_hz = s.bin_hz(peak_bin);
|
||||
assert!((peak_hz - freq as f32).abs() < 100.0, "peak at {peak_hz} Hz, want {freq}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn empty_input_yields_no_frames() {
|
||||
let s = analyze(&[], 48_000, 512, 256);
|
||||
assert!(s.frames.is_empty());
|
||||
}
|
||||
}
|
||||
+139
@@ -0,0 +1,139 @@
|
||||
//! A minimal canonical-PCM WAV reader (the analysis-side counterpart to
|
||||
//! `audio::recorder::WavWriter`). Handles 16-bit integer PCM — the only format
|
||||
//! PeerSpeak writes — in mono or interleaved multi-channel, downmixing to mono
|
||||
//! `f32` in `[-1, 1]` for analysis. I/O lives here, at the edge; the rest of the
|
||||
//! `dsp` module is pure.
|
||||
|
||||
use std::fs;
|
||||
use std::path::Path;
|
||||
|
||||
/// A decoded WAV: mono `f32` samples plus the original sample rate.
|
||||
pub struct WavData {
|
||||
pub samples: Vec<f32>,
|
||||
pub sample_rate: u32,
|
||||
}
|
||||
|
||||
/// Reads a 16-bit PCM WAV file, downmixing any channels to mono `f32`.
|
||||
///
|
||||
/// Returns `Err` with a human-readable reason if the file is missing, truncated,
|
||||
/// not RIFF/WAVE, or not 16-bit PCM. The parser walks the chunk list rather than
|
||||
/// assuming a fixed 44-byte header, so files with extra chunks (`LIST`, `fact`,
|
||||
/// …) still read.
|
||||
pub fn read(path: &Path) -> Result<WavData, String> {
|
||||
let bytes = fs::read(path).map_err(|e| format!("cannot read {}: {e}", path.display()))?;
|
||||
if bytes.len() < 12 || &bytes[0..4] != b"RIFF" || &bytes[8..12] != b"WAVE" {
|
||||
return Err("not a RIFF/WAVE file".to_string());
|
||||
}
|
||||
|
||||
let mut channels = 0u16;
|
||||
let mut sample_rate = 0u32;
|
||||
let mut bits = 0u16;
|
||||
let mut data: Option<&[u8]> = None;
|
||||
|
||||
// Walk chunks starting after the 12-byte RIFF/WAVE header.
|
||||
let mut pos = 12usize;
|
||||
while pos + 8 <= bytes.len() {
|
||||
let id = &bytes[pos..pos + 4];
|
||||
let size = u32::from_le_bytes([bytes[pos + 4], bytes[pos + 5], bytes[pos + 6], bytes[pos + 7]]) as usize;
|
||||
let body_start = pos + 8;
|
||||
let body_end = (body_start + size).min(bytes.len());
|
||||
match id {
|
||||
b"fmt " if size >= 16 => {
|
||||
let fmt = &bytes[body_start..body_end];
|
||||
let audio_format = u16::from_le_bytes([fmt[0], fmt[1]]);
|
||||
channels = u16::from_le_bytes([fmt[2], fmt[3]]);
|
||||
sample_rate = u32::from_le_bytes([fmt[4], fmt[5], fmt[6], fmt[7]]);
|
||||
bits = u16::from_le_bytes([fmt[14], fmt[15]]);
|
||||
if audio_format != 1 {
|
||||
return Err(format!("unsupported WAV format tag {audio_format} (need PCM=1)"));
|
||||
}
|
||||
}
|
||||
b"data" => {
|
||||
data = Some(&bytes[body_start..body_end]);
|
||||
}
|
||||
_ => {}
|
||||
}
|
||||
// Chunks are word-aligned: an odd size is padded with one byte.
|
||||
pos = body_start + size + (size & 1);
|
||||
}
|
||||
|
||||
if bits != 16 {
|
||||
return Err(format!("only 16-bit PCM supported, got {bits}-bit"));
|
||||
}
|
||||
let channels = channels.max(1);
|
||||
let data = data.ok_or("no data chunk")?;
|
||||
|
||||
// Interleaved S16LE -> per-frame channel average -> mono f32.
|
||||
let frame_bytes = 2 * channels as usize;
|
||||
let mut samples = Vec::with_capacity(data.len() / frame_bytes.max(1));
|
||||
for frame in data.chunks_exact(frame_bytes) {
|
||||
let mut acc = 0i32;
|
||||
for ch in frame.chunks_exact(2) {
|
||||
acc += i16::from_le_bytes([ch[0], ch[1]]) as i32;
|
||||
}
|
||||
let avg = acc as f32 / channels as f32;
|
||||
samples.push(avg / 32768.0);
|
||||
}
|
||||
|
||||
Ok(WavData { samples, sample_rate })
|
||||
}
|
||||
|
||||
/// Writes mono `f32` samples (clamped to `[-1, 1]`) as a 16-bit PCM WAV. Used by
|
||||
/// the `specview gen` command to materialize conjured test signals.
|
||||
pub fn write(path: &Path, samples: &[f32], sample_rate: u32) -> Result<(), String> {
|
||||
let pcm: Vec<i16> = samples
|
||||
.iter()
|
||||
.map(|&s| (s.clamp(-1.0, 1.0) * 32767.0).round() as i16)
|
||||
.collect();
|
||||
let data_bytes = (pcm.len() * 2) as u32;
|
||||
let byte_rate = sample_rate * 2; // mono, 2 bytes/sample
|
||||
let mut out = Vec::with_capacity(44 + pcm.len() * 2);
|
||||
out.extend_from_slice(b"RIFF");
|
||||
out.extend_from_slice(&(36 + data_bytes).to_le_bytes());
|
||||
out.extend_from_slice(b"WAVE");
|
||||
out.extend_from_slice(b"fmt ");
|
||||
out.extend_from_slice(&16u32.to_le_bytes());
|
||||
out.extend_from_slice(&1u16.to_le_bytes()); // PCM
|
||||
out.extend_from_slice(&1u16.to_le_bytes()); // mono
|
||||
out.extend_from_slice(&sample_rate.to_le_bytes());
|
||||
out.extend_from_slice(&byte_rate.to_le_bytes());
|
||||
out.extend_from_slice(&2u16.to_le_bytes()); // block align
|
||||
out.extend_from_slice(&16u16.to_le_bytes()); // bits
|
||||
out.extend_from_slice(b"data");
|
||||
out.extend_from_slice(&data_bytes.to_le_bytes());
|
||||
for s in pcm {
|
||||
out.extend_from_slice(&s.to_le_bytes());
|
||||
}
|
||||
fs::write(path, out).map_err(|e| format!("cannot write {}: {e}", path.display()))
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn write_then_read_round_trips() {
|
||||
let dir = std::env::temp_dir();
|
||||
let path = dir.join(format!("dsp-wav-rt-{}.wav", std::process::id()));
|
||||
let sig: Vec<f32> = (0..1000).map(|i| (i as f32 * 0.01).sin() * 0.5).collect();
|
||||
write(&path, &sig, 48_000).unwrap();
|
||||
|
||||
let back = read(&path).unwrap();
|
||||
assert_eq!(back.sample_rate, 48_000);
|
||||
assert_eq!(back.samples.len(), sig.len());
|
||||
// 16-bit quantization tolerance.
|
||||
for (a, b) in sig.iter().zip(&back.samples) {
|
||||
assert!((a - b).abs() < 1e-3, "{a} vs {b}");
|
||||
}
|
||||
let _ = fs::remove_file(&path);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn rejects_non_wave() {
|
||||
let dir = std::env::temp_dir();
|
||||
let path = dir.join(format!("dsp-wav-bad-{}.bin", std::process::id()));
|
||||
fs::write(&path, b"not a wav at all").unwrap();
|
||||
assert!(read(&path).is_err());
|
||||
let _ = fs::remove_file(&path);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,55 @@
|
||||
//! Analysis windows for the STFT.
|
||||
//!
|
||||
//! A raw rectangular frame leaks spectral energy across bins (the abrupt edges
|
||||
//! look like discontinuities to the FFT). A Hann window tapers each frame to
|
||||
//! zero at its edges, trading a little main-lobe width for much lower side-lobe
|
||||
//! leakage — the standard choice for a spectrogram.
|
||||
|
||||
use std::f64::consts::PI;
|
||||
|
||||
/// Returns a length-`n` periodic Hann window, `w[i] = 0.5·(1 - cos(2πi/n))`.
|
||||
///
|
||||
/// The *periodic* form (denominator `n`, not `n-1`) is used because STFT frames
|
||||
/// tile the signal; it gives perfect overlap-add reconstruction at 50% hop.
|
||||
pub fn hann(n: usize) -> Vec<f32> {
|
||||
if n <= 1 {
|
||||
return vec![1.0; n];
|
||||
}
|
||||
(0..n)
|
||||
.map(|i| (0.5 - 0.5 * (2.0 * PI * i as f64 / n as f64).cos()) as f32)
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Applies `window` to `frame` element-wise into a new buffer. Lengths must match.
|
||||
pub fn apply(frame: &[f32], window: &[f32]) -> Vec<f32> {
|
||||
debug_assert_eq!(frame.len(), window.len());
|
||||
frame.iter().zip(window).map(|(&s, &w)| s * w).collect()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn hann_endpoints_are_zero_and_center_is_one() {
|
||||
let w = hann(8);
|
||||
assert!(w[0].abs() < 1e-6, "first sample tapers to ~0, got {}", w[0]);
|
||||
// Periodic Hann peaks at the midpoint n/2.
|
||||
assert!((w[4] - 1.0).abs() < 1e-6, "center peaks at 1, got {}", w[4]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hann_is_symmetric_about_center() {
|
||||
let w = hann(16);
|
||||
// Periodic window is symmetric across indices 1..n-1.
|
||||
for i in 1..8 {
|
||||
assert!((w[i] - w[16 - i]).abs() < 1e-6);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn apply_scales_samples() {
|
||||
let out = apply(&[2.0, 2.0], &[0.5, 0.25]);
|
||||
assert_eq!(out, vec![1.0, 0.5]);
|
||||
}
|
||||
}
|
||||
@@ -1,5 +1,6 @@
|
||||
pub mod audio;
|
||||
pub mod codec;
|
||||
pub mod dsp;
|
||||
pub mod network;
|
||||
pub mod core;
|
||||
pub mod app;
|
||||
|
||||
Reference in New Issue
Block a user