mirror of
https://github.com/ruvnet/RuView
synced 2026-08-02 19:11:46 +00:00
perf(beyond-sota): ADR-154 M2 — FFT planner hoist (1.84x, bit-identical) + 3 honest perf nulls + boundary tests (#1055)
* perf(signal): hoist FFT planner across subcarriers (ADR-154 §7.4 #20) compute_multi_subcarrier_spectrogram called compute_spectrogram once per subcarrier, and each call built a fresh FftPlanner + re-planned the same length-window_size FFT. Hoist the plan + window out of the per-subcarrier loop via a new compute_spectrogram_with_plan core that takes a pre-planned Arc<dyn Fft> and pre-built window. compute_spectrogram delegates to it (unchanged behaviour); the multi-subcarrier path plans once and reuses. MEASURED-HOT (dsp_perf_bench, this box): at 56 subcarriers, window 128, fresh-planner-per-subcarrier 467.88 µs -> hoisted-plan 254.75 µs = 1.84x; window 256: 627.27 µs -> 448.39 µs = 1.40x. Plan-forward cost alone is ~1.86 µs (w128), x56 subcarriers ~= the removed delta. Output is bit-identical: multi_subcarrier_hoisted_plan_bit_identical compares f64::to_bits of every spectrogram value + freq/time resolution against the per-call fresh-planner path across all 4 window functions x {power,magnitude} on a 56-subcarrier matrix. The numeric STFT body is the old loop verbatim; only plan/window construction is lifted. Co-Authored-By: claude-flow <ruv@ruv.net> * test(signal): boundary/tolerance tests for ADR-154 §7.4 #14 #16 #19 Three "+ test" backlog gaps closed — pure additions, no behaviour change (phase_align refactor is internal: estimate_phase_offsets still returns the identical offset vector; a counted core is split out only to observe the iteration count). #14 cir.rs fft_operator — fft_operator_within_tolerance_of_dense_canonical56: the opt-in FFT Φ/Φᴴ path changes the witness hash, so pin it numerically CLOSE to the dense path (not silently divergent). Asserts the full Cir output (every tap within 1e-2·dominant, dominant idx/ratio, active_tap_count, ranging_valid, rms_delay_spread) on the production canonical-56 config across τ ∈ {20,50,90} ns. Extends the existing HT20/single-τ test. #16 phase_align.rs — refinement_terminates_at_iteration_cap_when_not_converging: forces non-convergence (tolerance=0.0, unreachable) and asserts the loop runs exactly max_iterations then returns — proving the cap, not convergence, bounds the loop (no infinite spin). Companion refinement_converges_before_cap_on_easy_input proves the cap is an upper bound, not the only exit. #19 csi_ratio.rs — ratio_finite_at_and_below_1e_12_epsilon: the module implements the CSI ratio as the conjugate product H_i·conj(H_j) (no division), so it is finite even at/below the 1e-12 magnitude boundary a naive H_i/H_j division would need an epsilon to guard. Pins finiteness + bit-exact conjugate product at the boundary (zero target → zero, never inf/NaN), through the amplitude/phase extraction. cargo test -p wifi-densepose-signal --no-default-features --lib: 447 passed, 0 failed; --features cir --lib: 447 passed, 0 failed. Co-Authored-By: claude-flow <ruv@ruv.net> * docs(adr-154): record Milestone-2 P2-perf verdicts + boundary tests (§7.4) §7.4: #20 MEASURED-HOT (1.40–1.84× spectrogram FFT-plan hoist, bit-identical); #5/#6/#7 MEASURED-NULL (benched, not hot, left as-is — sub-µs / stack-only / alloc-once); #8 MEASUREMENT-ONLY (per-call 56×56 eigh cost; eigenvalue/BLAS backend un-buildable on this Windows host, number deferred to a BLAS box, NOT fabricated; also corrects the finding — extract_perturbation reuses cached modes, the recompute is in estimate_occupancy). #14/#16/#19 RESOLVED (tolerance / convergence-cap / epsilon-boundary tests). Updated §7.4 intro + Horizon-ledger (deferred count 41→36). CHANGELOG [Unreleased] entry added. Co-Authored-By: claude-flow <ruv@ruv.net> * bench(signal): committed P2 bench-first benches (ADR-154 §7.4 #5/#6/#7/#8/#20) New dsp_perf_bench.rs backs every Milestone-2 perf verdict with a committed criterion bench — no speedup claimed without a before/after number here, and a benched NULL is the proof a micro-opt was unnecessary (the §5.x "already amortized" pattern). Registered in Cargo.toml [[bench]]. MEASURED (this box, criterion medians): #20 spectrogram_multi_subcarrier (fresh vs hoisted plan): MEASURED-HOT — 467.88→254.75 µs (1.84x) @ sc56/w128; 627.27→448.39 µs (1.40x) @ sc56/w256. Optimized in the prior commit. #5 multistatic_attention/weights: MEASURED-NULL — 181 ns (2 nodes) .. 848 ns (8 nodes); sub-µs, no hot-path alloc — left as-is. #6 tomography_reconstruct/solve: MEASURED-NULL — 47.5 µs (16 links) / 60.4 µs (32 links) for a full 50-iter ISTA solve; the 2 per-solve voxel buffers (~4 KB) are negligible vs O(iters·links·voxels) compute, and reconstruct(&self) reuses them across iterations already — left as-is. #7 pose_kalman_update/cycles: MEASURED-NULL — 150 ns (17 kpts) / 2.82 µs (170); the Kalman "gain matrices" are fixed-size STACK arrays ([[f32;3];6]), zero heap — nothing to reuse — left as-is. #8 field_model_occupancy (eigenvalue feature): MEASUREMENT-ONLY — quantifies the per-call n×n eigendecomposition cost; incremental SVD is a sized future project, not attempted (number recorded in ADR-154 §7.4). Reproduce: cargo bench -p wifi-densepose-signal --no-default-features --bench dsp_perf_bench cargo bench -p wifi-densepose-signal --bench dsp_perf_bench # adds #8 Cargo.lock: dev-dep (criterion/clap) graph + crate version bumps from the build; no runtime-dependency change. Co-Authored-By: claude-flow <ruv@ruv.net>
This commit is contained in:
rUv
Generated
+3
-3
@@ -10835,7 +10835,7 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "wifi-densepose-cli"
|
||||
version = "0.3.0"
|
||||
version = "0.3.1"
|
||||
dependencies = [
|
||||
"anyhow",
|
||||
"assert_cmd",
|
||||
@@ -11067,7 +11067,7 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "wifi-densepose-sensing-server"
|
||||
version = "0.3.2"
|
||||
version = "0.3.3"
|
||||
dependencies = [
|
||||
"axum",
|
||||
"chrono",
|
||||
@@ -11101,7 +11101,7 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "wifi-densepose-signal"
|
||||
version = "0.3.3"
|
||||
version = "0.3.4"
|
||||
dependencies = [
|
||||
"chrono",
|
||||
"criterion",
|
||||
|
||||
@@ -71,6 +71,12 @@ harness = false
|
||||
name = "features_bench"
|
||||
harness = false
|
||||
|
||||
## ADR-154 Milestone-2: P2 "bench-first" perf items (§7.4 #5/#6/#7/#8/#20).
|
||||
## #8 (field_model eigendecompose) is measured only under the eigenvalue feature.
|
||||
[[bench]]
|
||||
name = "dsp_perf_bench"
|
||||
harness = false
|
||||
|
||||
## ADR-134: CIR estimator throughput benchmarks
|
||||
[[bench]]
|
||||
name = "cir_bench"
|
||||
|
||||
@@ -0,0 +1,353 @@
|
||||
//! ADR-154 Milestone-2 perf benchmarks (§7.4 P2 "bench-first" items).
|
||||
//!
|
||||
//! PROOF discipline (ADR-154 §0): every P2 item is **benched before touched**.
|
||||
//! A micro-opt is landed only if the bench proves the path hot; otherwise the
|
||||
//! committed bench *is* the result — a MEASURED-NULL that proves the rewrite was
|
||||
//! unnecessary (exactly the §5.x "already amortized" pattern). No speedup is
|
||||
//! claimed without a before/after number from here.
|
||||
//!
|
||||
//! Reproduce (compile-only):
|
||||
//! cargo bench -p wifi-densepose-signal --no-default-features \
|
||||
//! --bench dsp_perf_bench --no-run
|
||||
//!
|
||||
//! Reproduce (full run, writes target/criterion/ HTML):
|
||||
//! cargo bench -p wifi-densepose-signal --no-default-features --bench dsp_perf_bench
|
||||
//!
|
||||
//! Groups:
|
||||
//! * `multistatic_attention` (#5) — `node_attention_weights` at 2..8 nodes ×
|
||||
//! 56 subcarriers. Re-derives consensus/softmax each call; no scratch to
|
||||
//! reuse → expected MEASURED-NULL.
|
||||
//! * `tomography_reconstruct` (#6) — full ISTA solve. The two voxel buffers are
|
||||
//! allocated once per `reconstruct()` (then `.fill`-reused across
|
||||
//! iterations), so the per-solve alloc is 2×n_voxels vs an
|
||||
//! O(iters·links·voxels) compute → expected MEASURED-NULL.
|
||||
//! * `pose_kalman_update` (#7) — Kalman predict+update loop. The "gain
|
||||
//! matrices" are fixed-size **stack** arrays (`[[f32;3];6]`), not heap —
|
||||
//! nothing to reuse → expected MEASURED-NULL.
|
||||
//! * `spectrogram_multi_subcarrier` (#20) — `compute_multi_subcarrier_spectrogram`:
|
||||
//! fresh-planner-per-subcarrier (BEFORE) vs hoisted-plan (AFTER, shipped).
|
||||
//! The per-subcarrier FFT re-plan is the likely real win.
|
||||
//! * `field_model_occupancy` (#8, `eigenvalue` only) — per-call n×n
|
||||
//! eigendecomposition in `estimate_occupancy`. MEASUREMENT-ONLY: quantifies
|
||||
//! the recompute cost; incremental SVD is a sized future project, not a
|
||||
//! micro-fix.
|
||||
|
||||
use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion, Throughput};
|
||||
use ndarray::Array2;
|
||||
use rustfft::FftPlanner;
|
||||
use std::f64::consts::PI;
|
||||
use std::time::Duration;
|
||||
|
||||
use wifi_densepose_signal::ruvsense::multistatic::node_attention_weights;
|
||||
use wifi_densepose_signal::ruvsense::pose_tracker::KeypointState;
|
||||
use wifi_densepose_signal::ruvsense::tomography::{
|
||||
LinkGeometry, Position3D, RfTomographer, TomographyConfig,
|
||||
};
|
||||
use wifi_densepose_signal::spectrogram::{
|
||||
compute_multi_subcarrier_spectrogram, compute_spectrogram, Spectrogram, SpectrogramConfig,
|
||||
WindowFunction,
|
||||
};
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// #5 multistatic node_attention_weights
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
fn make_node_amplitudes(n_nodes: usize, n_sub: usize) -> Vec<Vec<f32>> {
|
||||
(0..n_nodes)
|
||||
.map(|n| {
|
||||
(0..n_sub)
|
||||
.map(|s| {
|
||||
let phase = (n as f32 * 0.31 + s as f32 * 0.07) % std::f32::consts::TAU;
|
||||
0.5 + 0.4 * phase.sin()
|
||||
})
|
||||
.collect()
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn bench_multistatic_attention(c: &mut Criterion) {
|
||||
let mut group = c.benchmark_group("multistatic_attention");
|
||||
group.measurement_time(Duration::from_secs(3));
|
||||
let n_sub = 56; // canonical-56 grid
|
||||
|
||||
for &n_nodes in &[2usize, 4, 8] {
|
||||
let owned = make_node_amplitudes(n_nodes, n_sub);
|
||||
let refs: Vec<&[f32]> = owned.iter().map(|v| v.as_slice()).collect();
|
||||
group.throughput(Throughput::Elements(1));
|
||||
group.bench_with_input(
|
||||
BenchmarkId::new("weights", n_nodes),
|
||||
&refs,
|
||||
|b, amplitudes| {
|
||||
b.iter(|| black_box(node_attention_weights(black_box(amplitudes), 1.0)));
|
||||
},
|
||||
);
|
||||
}
|
||||
group.finish();
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// #6 tomography reconstruct (ISTA L1)
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
fn make_tomographer(n_links: usize) -> (RfTomographer, Vec<f64>) {
|
||||
// A modest 8x8x4 grid (256 voxels), n_links TX/RX pairs around the box.
|
||||
let config = TomographyConfig {
|
||||
nx: 8,
|
||||
ny: 8,
|
||||
nz: 4,
|
||||
bounds: [0.0, 0.0, 0.0, 4.0, 4.0, 2.0],
|
||||
lambda: 0.01,
|
||||
max_iterations: 50,
|
||||
tolerance: 1e-6,
|
||||
min_links: 8,
|
||||
};
|
||||
let mut links = Vec::with_capacity(n_links);
|
||||
for i in 0..n_links {
|
||||
let t = i as f64 / n_links as f64;
|
||||
links.push(LinkGeometry {
|
||||
tx: Position3D {
|
||||
x: 4.0 * (t * PI).cos().abs(),
|
||||
y: 0.0,
|
||||
z: 1.0,
|
||||
},
|
||||
rx: Position3D {
|
||||
x: 4.0 * (t * PI).sin().abs(),
|
||||
y: 4.0,
|
||||
z: 1.0,
|
||||
},
|
||||
link_id: i,
|
||||
});
|
||||
}
|
||||
let tomo = RfTomographer::new(config, &links).unwrap();
|
||||
// Deterministic attenuations (one occupied region in the middle).
|
||||
let attenuations: Vec<f64> = (0..n_links)
|
||||
.map(|i| 0.1 + 0.05 * ((i as f64 * 0.3).sin()))
|
||||
.collect();
|
||||
(tomo, attenuations)
|
||||
}
|
||||
|
||||
fn bench_tomography_reconstruct(c: &mut Criterion) {
|
||||
let mut group = c.benchmark_group("tomography_reconstruct");
|
||||
group.measurement_time(Duration::from_secs(4));
|
||||
|
||||
for &n_links in &[16usize, 32] {
|
||||
let (tomo, atten) = make_tomographer(n_links);
|
||||
group.throughput(Throughput::Elements(1));
|
||||
group.bench_with_input(
|
||||
BenchmarkId::new("solve", n_links),
|
||||
&(tomo, atten),
|
||||
|b, (tomo, atten)| {
|
||||
b.iter(|| black_box(tomo.reconstruct(black_box(atten)).unwrap().occupied_count));
|
||||
},
|
||||
);
|
||||
}
|
||||
group.finish();
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// #7 pose tracker Kalman update loop
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
fn bench_pose_kalman_update(c: &mut Criterion) {
|
||||
let mut group = c.benchmark_group("pose_kalman_update");
|
||||
group.measurement_time(Duration::from_secs(3));
|
||||
|
||||
// 17 keypoints (COCO-17), N predict+update cycles — a realistic frame batch.
|
||||
for &n_updates in &[17usize, 170] {
|
||||
group.throughput(Throughput::Elements(n_updates as u64));
|
||||
group.bench_with_input(BenchmarkId::new("cycles", n_updates), &n_updates, |b, &n| {
|
||||
b.iter(|| {
|
||||
let mut acc = 0.0_f32;
|
||||
for k in 0..n {
|
||||
let mut state = KeypointState::new(
|
||||
(k as f32 * 0.1).sin(),
|
||||
(k as f32 * 0.2).cos(),
|
||||
1.0 + (k as f32 * 0.05),
|
||||
);
|
||||
state.predict(0.05, 0.5);
|
||||
let meas = [
|
||||
(k as f32 * 0.1).sin() + 0.01,
|
||||
(k as f32 * 0.2).cos() - 0.01,
|
||||
1.0 + (k as f32 * 0.05),
|
||||
];
|
||||
state.update(&meas, 0.1, 1.0);
|
||||
acc += state.state[0];
|
||||
}
|
||||
black_box(acc)
|
||||
});
|
||||
});
|
||||
}
|
||||
group.finish();
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// #20 multi-subcarrier spectrogram: fresh-planner vs hoisted plan
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
fn make_csi_temporal(n_samples: usize, n_sc: usize) -> Array2<f64> {
|
||||
Array2::from_shape_fn((n_samples, n_sc), |(t, sc)| {
|
||||
let freq = 0.7 + sc as f64 * 0.13;
|
||||
(2.0 * PI * freq * t as f64 / 100.0).sin()
|
||||
+ 0.3 * (2.0 * PI * (freq * 2.1) * t as f64 / 100.0).cos()
|
||||
})
|
||||
}
|
||||
|
||||
/// BEFORE: re-plan the FFT inside `compute_spectrogram` for every subcarrier.
|
||||
/// Faithful transcription of the pre-ADR-154-M2 `compute_multi_subcarrier_spectrogram`.
|
||||
fn multi_fresh_planner(
|
||||
csi: &Array2<f64>,
|
||||
sample_rate: f64,
|
||||
config: &SpectrogramConfig,
|
||||
) -> Vec<Spectrogram> {
|
||||
let (_, n_sc) = csi.dim();
|
||||
(0..n_sc)
|
||||
.map(|sc| {
|
||||
let col: Vec<f64> = csi.column(sc).to_vec();
|
||||
// compute_spectrogram builds a fresh FftPlanner on every call.
|
||||
compute_spectrogram(&col, sample_rate, config).unwrap()
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn bench_spectrogram_multi_subcarrier(c: &mut Criterion) {
|
||||
let mut group = c.benchmark_group("spectrogram_multi_subcarrier");
|
||||
group.measurement_time(Duration::from_secs(5));
|
||||
let sample_rate = 100.0;
|
||||
|
||||
// Realistic: 600 temporal samples (~6 s @ 100 Hz) across 56 subcarriers,
|
||||
// window 128. n_sc re-plans removed by the hoist.
|
||||
for &(n_samples, n_sc, window) in &[(600usize, 56usize, 128usize), (600, 56, 256)] {
|
||||
let csi = make_csi_temporal(n_samples, n_sc);
|
||||
let config = SpectrogramConfig {
|
||||
window_size: window,
|
||||
hop_size: 64,
|
||||
window_fn: WindowFunction::Hann,
|
||||
power: true,
|
||||
};
|
||||
group.throughput(Throughput::Elements(n_sc as u64));
|
||||
|
||||
// BEFORE: fresh planner per subcarrier.
|
||||
group.bench_with_input(
|
||||
BenchmarkId::new("fresh_planner", format!("sc{n_sc}_w{window}")),
|
||||
&config,
|
||||
|b, cfg| {
|
||||
b.iter(|| black_box(multi_fresh_planner(black_box(&csi), sample_rate, cfg).len()));
|
||||
},
|
||||
);
|
||||
|
||||
// AFTER: hoisted plan (the shipped `compute_multi_subcarrier_spectrogram`).
|
||||
group.bench_with_input(
|
||||
BenchmarkId::new("hoisted_plan", format!("sc{n_sc}_w{window}")),
|
||||
&config,
|
||||
|b, cfg| {
|
||||
b.iter(|| {
|
||||
black_box(
|
||||
compute_multi_subcarrier_spectrogram(black_box(&csi), sample_rate, cfg)
|
||||
.unwrap()
|
||||
.len(),
|
||||
)
|
||||
});
|
||||
},
|
||||
);
|
||||
}
|
||||
group.finish();
|
||||
}
|
||||
|
||||
// A standalone FftPlanner sanity micro-bench documenting the cost the hoist
|
||||
// removes: building+planning a length-N forward FFT once.
|
||||
fn bench_fft_plan_cost(c: &mut Criterion) {
|
||||
let mut group = c.benchmark_group("fft_plan_cost");
|
||||
group.measurement_time(Duration::from_secs(2));
|
||||
for &n in &[128usize, 256] {
|
||||
group.bench_with_input(BenchmarkId::new("plan_forward", n), &n, |b, &n| {
|
||||
b.iter(|| {
|
||||
let mut planner = FftPlanner::<f64>::new();
|
||||
black_box(planner.plan_fft_forward(black_box(n)))
|
||||
});
|
||||
});
|
||||
}
|
||||
group.finish();
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// #8 field_model SVD/eigendecomposition recompute (MEASUREMENT-ONLY)
|
||||
// ---------------------------------------------------------------------------
|
||||
// `estimate_occupancy` builds an n×n covariance and eigendecomposes it on every
|
||||
// call (BLAS, `eigenvalue` feature). This bench quantifies that per-call cost so
|
||||
// ADR-154 §7.4 #8 can record a number; incremental SVD is a sized future item,
|
||||
// NOT attempted here.
|
||||
#[cfg(feature = "eigenvalue")]
|
||||
mod eig {
|
||||
use super::*;
|
||||
use wifi_densepose_signal::ruvsense::field_model::{FieldModel, FieldModelConfig};
|
||||
|
||||
fn calibrated_model(n_sub: usize, n_links: usize) -> FieldModel {
|
||||
let config = FieldModelConfig {
|
||||
n_subcarriers: n_sub,
|
||||
n_links,
|
||||
n_modes: 3,
|
||||
min_calibration_frames: 20,
|
||||
baseline_expiry_s: 86_400.0,
|
||||
};
|
||||
let mut model = FieldModel::new(config).unwrap();
|
||||
// Feed deterministic calibration frames: [n_links][n_sub] per observation.
|
||||
for f in 0..30 {
|
||||
let obs: Vec<Vec<f64>> = (0..n_links)
|
||||
.map(|l| {
|
||||
(0..n_sub)
|
||||
.map(|s| {
|
||||
0.5 + 0.3
|
||||
* ((f as f64 * 0.1 + l as f64 * 0.2 + s as f64 * 0.05).sin())
|
||||
})
|
||||
.collect()
|
||||
})
|
||||
.collect();
|
||||
model.feed_calibration(&obs).unwrap();
|
||||
}
|
||||
model.finalize_calibration(0, 0).unwrap();
|
||||
model
|
||||
}
|
||||
|
||||
pub fn bench_field_model_occupancy(c: &mut Criterion) {
|
||||
let mut group = c.benchmark_group("field_model_occupancy");
|
||||
group.measurement_time(Duration::from_secs(4));
|
||||
let n_sub = 56;
|
||||
let model = calibrated_model(n_sub, 4);
|
||||
// Sliding window of recent frames (50 ~ 2.5 s @ 20 Hz).
|
||||
let frames: Vec<Vec<f64>> = (0..50)
|
||||
.map(|t| {
|
||||
(0..n_sub)
|
||||
.map(|s| 0.5 + 0.3 * ((t as f64 * 0.15 + s as f64 * 0.07).sin()))
|
||||
.collect()
|
||||
})
|
||||
.collect();
|
||||
group.throughput(Throughput::Elements(1));
|
||||
group.bench_function(BenchmarkId::new("eigh", n_sub), |b| {
|
||||
b.iter(|| black_box(model.estimate_occupancy(black_box(&frames))));
|
||||
});
|
||||
group.finish();
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(feature = "eigenvalue")]
|
||||
criterion_group!(
|
||||
benches,
|
||||
bench_multistatic_attention,
|
||||
bench_tomography_reconstruct,
|
||||
bench_pose_kalman_update,
|
||||
bench_spectrogram_multi_subcarrier,
|
||||
bench_fft_plan_cost,
|
||||
eig::bench_field_model_occupancy,
|
||||
);
|
||||
|
||||
#[cfg(not(feature = "eigenvalue"))]
|
||||
criterion_group!(
|
||||
benches,
|
||||
bench_multistatic_attention,
|
||||
bench_tomography_reconstruct,
|
||||
bench_pose_kalman_update,
|
||||
bench_spectrogram_multi_subcarrier,
|
||||
bench_fft_plan_cost,
|
||||
);
|
||||
|
||||
criterion_main!(benches);
|
||||
@@ -197,4 +197,61 @@ mod tests {
|
||||
Err(CsiRatioError::LengthMismatch { .. })
|
||||
));
|
||||
}
|
||||
|
||||
// ADR-154 §7.4 #19: the CSI *ratio model*. The classic ratio is
|
||||
// `H_i[k] / H_j[k]`, which blows up (±inf / NaN) when `H_j[k]` approaches
|
||||
// zero — the case a `1e-12` division-guard epsilon is meant to protect. This
|
||||
// module deliberately implements the ratio as the **conjugate product**
|
||||
// `H_i * conj(H_j)` (SpotFi/IndoTrack), which has *no division* and is
|
||||
// therefore finite even at and below the `1e-12` magnitude boundary. This
|
||||
// test pins that property: at the epsilon boundary the output is finite and
|
||||
// exactly the conjugate product (no silent NaN/inf from a hidden divide).
|
||||
#[test]
|
||||
fn ratio_finite_at_and_below_1e_12_epsilon() {
|
||||
let eps = 1e-12_f64;
|
||||
// Reference at unit magnitude; target swept across / under the epsilon
|
||||
// boundary a naive H_i/H_j division would need to guard.
|
||||
let h_ref = vec![
|
||||
Complex64::from_polar(1.0, 0.3),
|
||||
Complex64::from_polar(1.0, 0.3),
|
||||
Complex64::from_polar(1.0, 0.3),
|
||||
Complex64::from_polar(1.0, 0.3),
|
||||
];
|
||||
let h_target = vec![
|
||||
Complex64::new(eps, 0.0), // exactly at the epsilon
|
||||
Complex64::new(eps * 0.5, 0.0), // below the epsilon
|
||||
Complex64::new(0.0, eps), // imaginary axis, at epsilon
|
||||
Complex64::new(0.0, 0.0), // exact zero — div would be inf/NaN
|
||||
];
|
||||
|
||||
let ratio = conjugate_multiply(&h_ref, &h_target).unwrap();
|
||||
assert_eq!(ratio.len(), 4);
|
||||
for (k, r) in ratio.iter().enumerate() {
|
||||
assert!(
|
||||
r.re.is_finite() && r.im.is_finite(),
|
||||
"conjugate-multiply ratio must be finite at boundary k={k}: {r:?}"
|
||||
);
|
||||
}
|
||||
|
||||
// The near-zero / zero target collapses the product toward zero (the
|
||||
// physically correct "no measurable path" answer), never to inf/NaN.
|
||||
assert!(
|
||||
ratio[3].norm() == 0.0,
|
||||
"exact-zero target → zero product, got {}",
|
||||
ratio[3].norm()
|
||||
);
|
||||
// The at-epsilon entries equal the exact conjugate product (bit-exact).
|
||||
let expected0 = h_ref[0] * h_target[0].conj();
|
||||
assert_eq!(ratio[0].re.to_bits(), expected0.re.to_bits());
|
||||
assert_eq!(ratio[0].im.to_bits(), expected0.im.to_bits());
|
||||
|
||||
// The full pipeline (amplitude/phase extraction) is also finite here.
|
||||
let mut m = Array2::<Complex64>::zeros((1, 4));
|
||||
for (k, &v) in ratio.iter().enumerate() {
|
||||
m[[0, k]] = v;
|
||||
}
|
||||
let (amp, phase) = ratio_to_amplitude_phase(&m);
|
||||
assert!(amp.iter().all(|a| a.is_finite()));
|
||||
assert!(phase.iter().all(|p| p.is_finite()));
|
||||
}
|
||||
}
|
||||
|
||||
@@ -1458,6 +1458,79 @@ mod tests {
|
||||
}
|
||||
}
|
||||
|
||||
/// ADR-154 §7.4 #14: the `fft_operator` path *changes the witness hash*
|
||||
/// (documented in `CirConfig::fft_operator`), so it must be pinned as
|
||||
/// numerically **close** to the dense path — not silently divergent. The
|
||||
/// existing `fft_estimate_matches_dense_dominant_tap` covers HT20 / one tau;
|
||||
/// this test asserts the **full `Cir` output** (every tap + every scalar
|
||||
/// field) stays within a documented relative tolerance on the production
|
||||
/// **canonical-56** config across several realistic delays. A regression
|
||||
/// that lets the FFT path drift (wrong scaling, off-by-one Φ column, etc.)
|
||||
/// fails here instead of corrupting a downstream witness unnoticed.
|
||||
#[test]
|
||||
fn fft_operator_within_tolerance_of_dense_canonical56() {
|
||||
// Relative tolerances — documented, not silent. The FFT operator sums the
|
||||
// same Φ entries in a different order, so taps agree to ~float epsilon
|
||||
// scaled by the dominant-tap magnitude; ISTA can differ by a few last
|
||||
// bits over its trajectory, hence 1e-2 (same order as the existing test).
|
||||
const TAP_REL_TOL: f32 = 1e-2;
|
||||
const RATIO_ABS_TOL: f32 = 1e-2;
|
||||
const SPREAD_REL_TOL: f64 = 1e-2;
|
||||
|
||||
for &tau in &[20e-9_f64, 50e-9, 90e-9] {
|
||||
let dense_cfg = CirConfig::canonical56();
|
||||
let mut fft_cfg = CirConfig::canonical56();
|
||||
fft_cfg.fft_operator = true;
|
||||
|
||||
let frame = make_single_tap_frame(dense_cfg.num_subcarriers, tau);
|
||||
let dense = CirEstimator::new(dense_cfg).estimate(&frame).unwrap();
|
||||
let fast = CirEstimator::new(fft_cfg).estimate(&frame).unwrap();
|
||||
|
||||
assert_eq!(dense.taps.len(), fast.taps.len());
|
||||
|
||||
// Full tap vector close (relative to the dominant tap magnitude).
|
||||
let dom = dense.taps[dense.dominant_tap_idx].norm().max(1e-6);
|
||||
let mut max_tap_err = 0.0_f32;
|
||||
for (a, b) in dense.taps.iter().zip(&fast.taps) {
|
||||
max_tap_err = max_tap_err.max((a - b).norm());
|
||||
}
|
||||
assert!(
|
||||
max_tap_err <= TAP_REL_TOL * dom,
|
||||
"tau={tau:e}: FFT taps diverged from dense — max err {max_tap_err} > {TAP_REL_TOL} * {dom} (NOT numerically close)"
|
||||
);
|
||||
|
||||
// The dominant tap and the scalar summary fields must agree too —
|
||||
// these feed the witness, so a silent divergence here is the bug #14
|
||||
// guards against.
|
||||
assert_eq!(
|
||||
dense.dominant_tap_idx, fast.dominant_tap_idx,
|
||||
"tau={tau:e}: dominant tap index moved"
|
||||
);
|
||||
assert!(
|
||||
(dense.dominant_tap_ratio - fast.dominant_tap_ratio).abs() <= RATIO_ABS_TOL,
|
||||
"tau={tau:e}: dominant_tap_ratio drift {} vs {}",
|
||||
dense.dominant_tap_ratio,
|
||||
fast.dominant_tap_ratio
|
||||
);
|
||||
assert_eq!(
|
||||
dense.active_tap_count, fast.active_tap_count,
|
||||
"tau={tau:e}: active_tap_count changed"
|
||||
);
|
||||
assert_eq!(
|
||||
dense.ranging_valid, fast.ranging_valid,
|
||||
"tau={tau:e}: ranging_valid flipped"
|
||||
);
|
||||
let spread_ref = dense.rms_delay_spread_s.abs().max(1e-12);
|
||||
assert!(
|
||||
(dense.rms_delay_spread_s - fast.rms_delay_spread_s).abs()
|
||||
<= SPREAD_REL_TOL * spread_ref,
|
||||
"tau={tau:e}: rms_delay_spread drift {} vs {}",
|
||||
dense.rms_delay_spread_s,
|
||||
fast.rms_delay_spread_s
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// The default configs keep the FFT operator off — the dense, bit-exact
|
||||
/// witness path is the default (enabling FFT shifts float results).
|
||||
#[test]
|
||||
|
||||
@@ -201,12 +201,29 @@ fn find_static_subcarriers(
|
||||
|
||||
/// Estimate per-channel phase offsets using iterative Neumann-style refinement.
|
||||
///
|
||||
/// Channel 0 is the reference (offset = 0).
|
||||
/// Channel 0 is the reference (offset = 0). Thin wrapper that drops the
|
||||
/// iteration count; `estimate_phase_offsets_counted` is the instrumented core.
|
||||
fn estimate_phase_offsets(
|
||||
frames: &[CanonicalCsiFrame],
|
||||
static_indices: &[usize],
|
||||
config: &PhaseAlignConfig,
|
||||
) -> std::result::Result<Vec<f32>, PhaseAlignError> {
|
||||
estimate_phase_offsets_counted(frames, static_indices, config).map(|(offsets, _iters)| offsets)
|
||||
}
|
||||
|
||||
/// Core of [`estimate_phase_offsets`], also returning the number of refinement
|
||||
/// iterations actually executed.
|
||||
///
|
||||
/// The returned count is bounded by `config.max_iterations` — that bound is the
|
||||
/// convergence cap that guarantees termination on inputs the damped Neumann
|
||||
/// update never drives below `config.tolerance` (ADR-154 §7.4 #16). The offset
|
||||
/// vector is identical to the public `estimate_phase_offsets` path; only the
|
||||
/// iteration count is surfaced (for the cap test).
|
||||
fn estimate_phase_offsets_counted(
|
||||
frames: &[CanonicalCsiFrame],
|
||||
static_indices: &[usize],
|
||||
config: &PhaseAlignConfig,
|
||||
) -> std::result::Result<(Vec<f32>, usize), PhaseAlignError> {
|
||||
let n_ch = frames.len();
|
||||
let mut offsets = vec![0.0_f32; n_ch];
|
||||
|
||||
@@ -220,7 +237,7 @@ fn estimate_phase_offsets(
|
||||
}
|
||||
|
||||
// Iterative refinement (Neumann-style)
|
||||
for _iter in 0..config.max_iterations {
|
||||
for iter in 0..config.max_iterations {
|
||||
let mut max_update = 0.0_f32;
|
||||
|
||||
for c in 1..n_ch {
|
||||
@@ -241,12 +258,13 @@ fn estimate_phase_offsets(
|
||||
}
|
||||
|
||||
if max_update < config.tolerance {
|
||||
return Ok(offsets);
|
||||
return Ok((offsets, iter + 1));
|
||||
}
|
||||
}
|
||||
|
||||
// Even if we do not converge tightly, return best estimate
|
||||
Ok(offsets)
|
||||
// Even if we do not converge tightly, return best estimate. The loop ran the
|
||||
// full cap — termination is guaranteed by `config.max_iterations`.
|
||||
Ok((offsets, config.max_iterations))
|
||||
}
|
||||
|
||||
/// Apply phase correction: subtract offset from each subcarrier phase.
|
||||
@@ -446,6 +464,73 @@ mod tests {
|
||||
assert_eq!(cfg.min_static_subcarriers, 5);
|
||||
}
|
||||
|
||||
// ADR-154 §7.4 #16: the iterative LO-offset refinement must TERMINATE at the
|
||||
// `max_iterations` cap on a non-converging input — no unbounded loop.
|
||||
//
|
||||
// We force non-convergence by setting `tolerance` to an unreachable value
|
||||
// (the damped Neumann update on bounded phase residuals can never drive
|
||||
// `max_update` below 0.0), so the `max_update < tolerance` early-exit is
|
||||
// never taken. The instrumented core must then run *exactly*
|
||||
// `max_iterations` and return — proving the cap, not convergence, is what
|
||||
// bounds the loop.
|
||||
#[test]
|
||||
fn refinement_terminates_at_iteration_cap_when_not_converging() {
|
||||
let n_sub = 56;
|
||||
let max_iterations = 7;
|
||||
let config = PhaseAlignConfig {
|
||||
max_iterations,
|
||||
// Unreachable tolerance: `max_update` is always ≥ 0, never < 0.0,
|
||||
// so the convergence branch can never fire.
|
||||
tolerance: 0.0,
|
||||
static_fraction: 0.3,
|
||||
min_static_subcarriers: 5,
|
||||
};
|
||||
// Two channels with a real, persistent offset so each iteration keeps
|
||||
// producing a non-zero update.
|
||||
let f0 = make_frame_with_phase(n_sub, 0.0, 0.0);
|
||||
let f1 = make_frame_with_phase(n_sub, 0.0, 1.3);
|
||||
let frames = vec![f0, f1];
|
||||
let static_indices = find_static_subcarriers(&frames, &config).unwrap();
|
||||
|
||||
let (offsets, iters) =
|
||||
estimate_phase_offsets_counted(&frames, &static_indices, &config).unwrap();
|
||||
|
||||
// The cap, not convergence, terminated the loop.
|
||||
assert_eq!(
|
||||
iters, max_iterations,
|
||||
"expected the loop to run the full cap ({max_iterations}), got {iters}"
|
||||
);
|
||||
// It still returns a finite best-estimate offset vector.
|
||||
assert_eq!(offsets.len(), 2);
|
||||
assert!(offsets.iter().all(|o| o.is_finite()));
|
||||
// Reference channel offset stays 0.
|
||||
assert_eq!(offsets[0], 0.0);
|
||||
}
|
||||
|
||||
// Convergent companion: a near-identical input converges *before* the cap,
|
||||
// so the cap is an upper bound, not the only exit.
|
||||
#[test]
|
||||
fn refinement_converges_before_cap_on_easy_input() {
|
||||
let n_sub = 56;
|
||||
let config = PhaseAlignConfig {
|
||||
max_iterations: 50,
|
||||
tolerance: 1e-2, // loose: a tiny offset converges in a few iters
|
||||
static_fraction: 0.3,
|
||||
min_static_subcarriers: 5,
|
||||
};
|
||||
let f0 = make_frame_with_phase(n_sub, 0.0, 0.0);
|
||||
let f1 = make_frame_with_phase(n_sub, 0.0, 0.02);
|
||||
let frames = vec![f0, f1];
|
||||
let static_indices = find_static_subcarriers(&frames, &config).unwrap();
|
||||
let (_offsets, iters) =
|
||||
estimate_phase_offsets_counted(&frames, &static_indices, &config).unwrap();
|
||||
assert!(
|
||||
iters < config.max_iterations,
|
||||
"easy input should converge before the cap, ran {iters}/{}",
|
||||
config.max_iterations
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn phase_correction_preserves_amplitude() {
|
||||
let mut aligner = PhaseAligner::new(2);
|
||||
|
||||
@@ -9,9 +9,10 @@
|
||||
|
||||
use ndarray::Array2;
|
||||
use num_complex::Complex64;
|
||||
use rustfft::FftPlanner;
|
||||
use rustfft::{Fft, FftPlanner};
|
||||
use ruvector_attn_mincut::attn_mincut;
|
||||
use std::f64::consts::PI;
|
||||
use std::sync::Arc;
|
||||
|
||||
/// Configuration for spectrogram generation.
|
||||
#[derive(Debug, Clone)]
|
||||
@@ -87,12 +88,40 @@ pub fn compute_spectrogram(
|
||||
return Err(SpectrogramError::InvalidWindowSize);
|
||||
}
|
||||
|
||||
let n_frames = (signal.len() - config.window_size) / config.hop_size + 1;
|
||||
let n_freq = config.window_size / 2 + 1;
|
||||
let window = make_window(config.window_fn, config.window_size);
|
||||
|
||||
let mut planner = FftPlanner::new();
|
||||
let fft = planner.plan_fft_forward(config.window_size);
|
||||
let window = make_window(config.window_fn, config.window_size);
|
||||
Ok(compute_spectrogram_with_plan(
|
||||
signal,
|
||||
sample_rate,
|
||||
config,
|
||||
&fft,
|
||||
&window,
|
||||
))
|
||||
}
|
||||
|
||||
/// STFT core that runs against a **pre-planned** FFT and pre-built window.
|
||||
///
|
||||
/// ADR-154 §7.4 #20: `compute_spectrogram` re-plans the FFT on every call, so
|
||||
/// `compute_multi_subcarrier_spectrogram` (which calls it once per subcarrier)
|
||||
/// re-planned the same length-`window_size` FFT for *every* subcarrier. This
|
||||
/// helper hoists the plan + window out of the per-subcarrier loop. The numeric
|
||||
/// body is byte-for-byte the old loop — only the plan/window construction is
|
||||
/// lifted — so the output is **bit-identical** to the per-call path (asserted by
|
||||
/// `multi_subcarrier_hoisted_plan_bit_identical`). Callers must pass a plan
|
||||
/// built for exactly `config.window_size` and a window of that length.
|
||||
fn compute_spectrogram_with_plan(
|
||||
signal: &[f64],
|
||||
sample_rate: f64,
|
||||
config: &SpectrogramConfig,
|
||||
fft: &Arc<dyn Fft<f64>>,
|
||||
window: &[f64],
|
||||
) -> Spectrogram {
|
||||
debug_assert_eq!(window.len(), config.window_size, "window/plan size mismatch");
|
||||
debug_assert_eq!(fft.len(), config.window_size, "FFT/window size mismatch");
|
||||
|
||||
let n_frames = (signal.len() - config.window_size) / config.hop_size + 1;
|
||||
let n_freq = config.window_size / 2 + 1;
|
||||
|
||||
let mut data = Array2::zeros((n_freq, n_frames));
|
||||
|
||||
@@ -116,13 +145,13 @@ pub fn compute_spectrogram(
|
||||
}
|
||||
}
|
||||
|
||||
Ok(Spectrogram {
|
||||
Spectrogram {
|
||||
data,
|
||||
n_freq,
|
||||
n_time: n_frames,
|
||||
freq_resolution: sample_rate / config.window_size as f64,
|
||||
time_resolution: config.hop_size as f64 / sample_rate,
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
/// Compute spectrogram for each subcarrier from a temporal CSI matrix.
|
||||
@@ -134,12 +163,40 @@ pub fn compute_multi_subcarrier_spectrogram(
|
||||
sample_rate: f64,
|
||||
config: &SpectrogramConfig,
|
||||
) -> Result<Vec<Spectrogram>, SpectrogramError> {
|
||||
let (_, n_sc) = csi_temporal.dim();
|
||||
let mut spectrograms = Vec::with_capacity(n_sc);
|
||||
let (n_samples, n_sc) = csi_temporal.dim();
|
||||
|
||||
// ADR-154 §7.4 #20: validate *once* (same checks `compute_spectrogram`
|
||||
// makes), then plan the FFT + build the window *once* and reuse them across
|
||||
// every subcarrier instead of re-planning per column. The window length is
|
||||
// identical for all subcarriers, so this is pure hoisting — output stays
|
||||
// bit-identical to the per-call path.
|
||||
if n_samples < config.window_size {
|
||||
return Err(SpectrogramError::SignalTooShort {
|
||||
signal_len: n_samples,
|
||||
window_size: config.window_size,
|
||||
});
|
||||
}
|
||||
if config.hop_size == 0 {
|
||||
return Err(SpectrogramError::InvalidHopSize);
|
||||
}
|
||||
if config.window_size == 0 {
|
||||
return Err(SpectrogramError::InvalidWindowSize);
|
||||
}
|
||||
|
||||
let mut planner = FftPlanner::new();
|
||||
let fft = planner.plan_fft_forward(config.window_size);
|
||||
let window = make_window(config.window_fn, config.window_size);
|
||||
|
||||
let mut spectrograms = Vec::with_capacity(n_sc);
|
||||
for sc in 0..n_sc {
|
||||
let col: Vec<f64> = csi_temporal.column(sc).to_vec();
|
||||
spectrograms.push(compute_spectrogram(&col, sample_rate, config)?);
|
||||
spectrograms.push(compute_spectrogram_with_plan(
|
||||
&col,
|
||||
sample_rate,
|
||||
config,
|
||||
&fft,
|
||||
&window,
|
||||
));
|
||||
}
|
||||
|
||||
Ok(spectrograms)
|
||||
@@ -372,6 +429,67 @@ mod tests {
|
||||
assert_eq!(spec.n_freq, 65);
|
||||
}
|
||||
}
|
||||
|
||||
// ADR-154 §7.4 #20: the FFT-planner hoist in
|
||||
// `compute_multi_subcarrier_spectrogram` must produce **bit-identical**
|
||||
// output to calling `compute_spectrogram` (fresh planner) per subcarrier.
|
||||
// We compare `f64::to_bits` of every spectrogram value across several
|
||||
// window functions and a realistic 56-subcarrier CSI matrix — the planner
|
||||
// change only reorders *when* the (identical) plan is built, never the math.
|
||||
#[test]
|
||||
fn multi_subcarrier_hoisted_plan_bit_identical() {
|
||||
let n_samples = 600;
|
||||
let n_sc = 56; // canonical-56 grid — the production subcarrier count
|
||||
let sample_rate = 100.0;
|
||||
let csi = Array2::from_shape_fn((n_samples, n_sc), |(t, sc)| {
|
||||
// Deterministic, non-trivial per-subcarrier content.
|
||||
let freq = 0.7 + sc as f64 * 0.13;
|
||||
(2.0 * PI * freq * t as f64 / sample_rate).sin()
|
||||
+ 0.3 * (2.0 * PI * (freq * 2.1) * t as f64 / sample_rate).cos()
|
||||
});
|
||||
|
||||
for window_fn in [
|
||||
WindowFunction::Hann,
|
||||
WindowFunction::Hamming,
|
||||
WindowFunction::Blackman,
|
||||
WindowFunction::Rectangular,
|
||||
] {
|
||||
for &power in &[true, false] {
|
||||
let config = SpectrogramConfig {
|
||||
window_size: 128,
|
||||
hop_size: 37, // non-divisor hop to exercise frame edges
|
||||
window_fn,
|
||||
power,
|
||||
};
|
||||
|
||||
// AFTER: hoisted-plan path.
|
||||
let hoisted =
|
||||
compute_multi_subcarrier_spectrogram(&csi, sample_rate, &config).unwrap();
|
||||
|
||||
// BEFORE: independent per-subcarrier fresh-planner path.
|
||||
let reference: Vec<Spectrogram> = (0..n_sc)
|
||||
.map(|sc| {
|
||||
let col: Vec<f64> = csi.column(sc).to_vec();
|
||||
compute_spectrogram(&col, sample_rate, &config).unwrap()
|
||||
})
|
||||
.collect();
|
||||
|
||||
assert_eq!(hoisted.len(), reference.len());
|
||||
for (sc, (h, r)) in hoisted.iter().zip(reference.iter()).enumerate() {
|
||||
assert_eq!(h.data.dim(), r.data.dim(), "dim sc={sc} {window_fn:?}");
|
||||
for (a, b) in h.data.iter().zip(r.data.iter()) {
|
||||
assert_eq!(
|
||||
a.to_bits(),
|
||||
b.to_bits(),
|
||||
"bit mismatch sc={sc} {window_fn:?} power={power}: {a} vs {b}"
|
||||
);
|
||||
}
|
||||
assert_eq!(h.freq_resolution.to_bits(), r.freq_resolution.to_bits());
|
||||
assert_eq!(h.time_resolution.to_bits(), r.time_resolution.to_bits());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
|
||||
Reference in New Issue
Block a user