Skip to main content

rill_core_model/
analysis.rs

1use crate::WdfElement;
2use parking_lot::RwLock;
3use rill_core::Transcendental;
4use std::sync::Arc;
5
6/// Analyze magnitude response of a WDF element chain
7///
8/// Returns a vector of (frequency, magnitude) pairs.
9pub fn frequency_response_magnitude<T: Transcendental>(
10    elements: &[Arc<RwLock<dyn WdfElement<T>>>],
11    frequencies: &[T],
12    _sample_rate: T,
13) -> Vec<(T, T)> {
14    let mut response = Vec::new();
15
16    for &freq in frequencies {
17        let two = T::from_f32(2.0);
18        let omega = two * T::PI * freq;
19        let mut mag = T::ONE;
20
21        for element in elements {
22            let r = element.read().port_resistance();
23            let denom_sq = T::ONE + omega * omega * r * r;
24            mag /= denom_sq.sqrt();
25        }
26
27        response.push((freq, mag));
28    }
29
30    response
31}
32
33/// Analyze total harmonic distortion (THD) of a WDF element
34pub fn analyze_distortion<T: Transcendental>(
35    element: &mut dyn WdfElement<T>,
36    frequency: T,
37    amplitude: T,
38    sample_rate: T,
39    num_cycles: usize,
40) -> T {
41    let num_samples = (sample_rate.to_f64() / frequency.to_f64() * num_cycles as f64) as usize;
42    let mut output = Vec::with_capacity(num_samples);
43
44    for i in 0..num_samples {
45        let two = T::from_f32(2.0);
46        let t = T::from_f32(i as f32) / sample_rate;
47        let sample = amplitude * (two * T::PI * frequency * t).sin();
48        let b = element.process_incident(sample);
49        output.push(b);
50        element.update_state();
51    }
52
53    let peak_output = output.iter().cloned().fold(T::ZERO, T::max);
54    let fundamental_amplitude = amplitude;
55
56    ((peak_output - fundamental_amplitude) / fundamental_amplitude).abs()
57}