lasprs 0.14.1

Library for Acoustic Signal Processing (Rust edition, with optional Python bindings via pyo3)
//! Computation of sound levels vs time over a measurement, using the
//! pre-capture data for the time weighting filter warm-up.
//!
//! The pre-capture dataset (the last 3 seconds before the recording
//! started) is used to let the SLM's time weighting filter settle, so
//! that the levels-vs-time output is valid from the start of the live
//! data (t = 0). When the warm-up time exceeds the pre-capture duration
//! (or when no pre-capture is present), the output starts later, at
//! `max(0, warmup - precap) / fs` seconds into the live data.
//!
//! This module is a thin wrapper: it feeds the measurement's time data
//! (pre-capture first, then live) to an [SLM] and returns the
//! accumulated [SLMResult]. All SLM logic (warm-up, negative-time
//! suppression, decimation) lives in the SLM itself.

use super::*;
use crate::filter::StandardFilterDescriptor;
use crate::slm::{SLM, SLMResult, SLMSettings};
use anyhow::{Result, anyhow};
use std::collections::BTreeMap;

impl SharedMeasurement {
    /// Compute sound levels vs time for a single channel of this
    /// measurement, on a background thread.
    ///
    /// The pre-capture data is used to warm up the time weighting
    /// filter, so that the levels-vs-time output is valid from the
    /// start of the live data (t = 0). The SLM statistics (Leq, Lmax,
    /// Lpk) are reset once the warm-up is over, so they only cover the
    /// valid period.
    ///
    /// # Arguments
    ///
    /// * `settings` - SLM settings (band descriptors, frequency and
    ///   time weighting, Lref, decimation)
    /// * `channel` - Input channel to compute levels for
    /// * `istart` - Start frame of the live data to include (defaults
    ///   to the beginning)
    /// * `istop` - Exclusive end frame of the live data (defaults to
    ///   the end)
    /// * `only_stats` - Only compute the statistics (Leq, Lmax, Lpk),
    ///   skip the levels-vs-time output
    pub fn SLM(
        &self,
        settings: SLMSettings,
        channel: usize,
        istart: Option<usize>,
        istop: Option<usize>,
        only_stats: bool,
    ) -> Result<SLMResult> {
        let meas = self.read();
        let nframes = meas.nframes();

        // The frame range and channel bounds are validated when the block
        // iterators are created below (ConvertedBlockIter / RawBlockIter).
        let istart = istart.unwrap_or(0);
        let istop = istop.unwrap_or(nframes);

        // Number of available pre-capture frames
        let total_precap_frames = meas.precapture_nframes() as i64;

        // Required SLM warmup time
        let warmup_time = settings.timeWeighting.warmupTime();

        // Ideal number of warmup frames (based on warmup time and samplerate)
        let ideal_warmpup_frames = (warmup_time.as_secs_f64() * *meas.samplerate()) as i64;

        let ideal_initial_frame = istart as i64 - ideal_warmpup_frames;
        let initial_frame = (ideal_initial_frame).max(-total_precap_frames);

        let mut slm = SLM::new(settings, Some(initial_frame));

        let mut lt_accum: Option<BTreeMap<_, _>> = None;
        let mut t_accum: Vec<Flt> = Vec::new();
        let mut last: Option<SLMResult> = None;

        // Create the live data iterator before possibly the pre-capture
        // iterator, for the error messages, that need to be reported first.
        let live_iter = meas.data_iter(
            Some(&[channel]),
            Some(initial_frame.max(0) as usize),
            Some(istop),
        )?;

        // Iterate over the pre-capture data first, then the live data. The
        // iterators borrow the measurement, so the read guard must stay
        // alive while iterating.`
        let precap_iter = if initial_frame < 0 {
            meas.precapture_iter(
                Some(&[channel]),
                Some((total_precap_frames + initial_frame) as usize),
                None,
            )?
        } else {
            None
        };

        for block in precap_iter.into_iter().flatten().chain(live_iter) {
            // Block should have only one column, as only one channel is selected
            debug_assert!(block.ncols() == 1, "Invalid number of columns");
            let td = block.column(0);
            let td = td.as_slice().expect("converted block must be contiguous");
            if let Some(result) = slm.run(td, !only_stats) {
                last = Some(result.clone());
                if let Some(lt) = &result.Lt {
                    let acc = lt_accum
                        .get_or_insert_with(|| lt.keys().map(|k| (*k, Vec::new())).collect());
                    for (desc, band) in lt {
                        acc.get_mut(desc)
                            .expect("accumulator exists for descriptor")
                            .extend_from_slice(band);
                    }
                }
                if let Some(t) = &result.t {
                    t_accum.extend(t.iter());
                }
            }
        }

        let mut result = last.ok_or_else(|| {
            anyhow!(
                "Not enough data for the SLM to produce valid output with the \
             chosen settings (warm-up never completed)."
            )
        })?;

        result.Lt = if only_stats { None } else { lt_accum };
        result.t = if t_accum.is_empty() {
            None
        } else {
            Some(t_accum)
        };
        Ok(result)
    }
}

#[cfg(test)]
mod test {
    use super::*;
    use crate::daq::RecordStatus;
    use crate::daq::{DaqChannel, DaqConfig, RecordSettings, Recording, StreamMgr, StreamType};
    use crate::filter::StandardFilterDescriptor;
    use crate::measurement::MeasurementType;
    use crate::siggen::{SiggenCommand, SourceDescriptor};
    use crate::slm::SLMSettingsBuilder;
    use crate::{FreqWeighting, TimeWeighting};
    use approx::assert_abs_diff_eq;
    use std::time::Duration;

    /// Record a sine tone via loopback and verify levels-vs-time output.
    /// Requires the `loopback-api` feature (enabled via `test_features`).
    #[cfg(feature = "loopback-api")]
    #[test]
    fn test_slm_on_measurement() -> anyhow::Result<()> {
        let tmpdir = tempfile::tempdir()?;
        let mut smgr = StreamMgr::new_with_devices();

        let in_cfg = DaqConfig {
            api: crate::daq::DaqApiDescriptor::Loopback,
            device_name: "Loopback".into(),
            inchannel_config: vec![
                DaqChannel::defaultAudio("ch0"),
                DaqChannel::defaultAudio("ch1"),
            ],
            outchannel_config: vec![],
            dtype: crate::daq::DataType::F32,
            sampleRateIndex: 0,
            framesPerBlockIndex: 0,
            monitorOutput: false,
        };
        let out_cfg = {
            let mut cfg = in_cfg.clone();
            cfg.inchannel_config = vec![];
            cfg.outchannel_config = vec![
                DaqChannel::defaultAudio("ch0"),
                DaqChannel::defaultAudio("ch1"),
            ];
            cfg
        };

        smgr.startStream(StreamType::Input, &in_cfg)?;
        smgr.startStream(StreamType::Output, &out_cfg)?;
        smgr.setSiggenSource(SourceDescriptor::Sine {
            frequency: (1000.).try_into().unwrap(),
        })?;
        smgr.siggenCommand(SiggenCommand::SetMuteAllChannels { mute: false })?;
        // Set amplitude to √2 → RMS = 1 → 94 dB SPL
        smgr.siggenCommand(SiggenCommand::SetAllGains {
            g: sq2.try_into().unwrap(),
        })?;

        let duration = Duration::from_secs(3);
        let rec_settings = RecordSettings::new(
            "slm_test",
            false,
            Some(duration),
            None,
            None,
            // Do not perform pre-capture
            Some(MeasurementType::NotSpecific {}),
            false,
            Some(tmpdir.path()),
            vec![],
        )?;

        let mut recording = Recording::new(rec_settings, &mut smgr)?;
        loop {
            use crate::daq::RecordStatus;

            if matches!(recording.status(), RecordStatus::Finished { .. }) {
                break;
            }
            std::thread::sleep(Duration::from_millis(10));
        }
        recording.stop();

        let m = recording.getMeasurement().unwrap();
        drop(recording);
        drop(smgr);

        let fs = *m.read().samplerate();
        let fsp = StrictlyPositive::new(fs)?;
        let desc = StandardFilterDescriptor::Overall().unwrap();
        let slm_settings = SLMSettingsBuilder::default()
            .fs(fsp)
            .timeWeighting(TimeWeighting::Fast {})
            .freqWeighting(FreqWeighting::Z)
            .filterDescriptors([desc])
            .Lref(StrictlyPositive::new(2e-5).unwrap())
            .build()
            .unwrap();

        // only_stats = true
        let result = m.SLM(slm_settings.clone(), 0, None, None, true).unwrap();

        // 1 kHz sine at √2 amplitude → RMS = 1 → Leq ≈ 93.98 dB SPL
        // Peak squared = 2 → Lpk ≈ 96.99 dB SPL
        // Measured: Lmax=93.9821, Lpk=96.9897, Leq=93.9794. The Lmax deficit
        // (~0.018 dB) is the Fast time-weighting asymptote, not measurement
        // noise. Epsilons leave ~10-30x headroom over measured deviations.
        assert_abs_diff_eq!(result.Lmax[&desc], 94.0, epsilon = 0.03);
        assert_abs_diff_eq!(result.Lpk[&desc], 96.99, epsilon = 0.001);
        assert_abs_diff_eq!(result.Leq[&desc], 93.98, epsilon = 0.002);
        assert!(
            result.Lt.is_none(),
            "Lt should be None when only_stats=true"
        );
        assert!(result.t.is_none());

        // only_stats = false: Lt should be present with time axis
        let result2 = m.SLM(slm_settings, 0, None, None, false).unwrap();
        assert_abs_diff_eq!(result2.Lmax[&desc], 94.0, epsilon = 0.03);
        assert_abs_diff_eq!(result2.Lpk[&desc], 96.99, epsilon = 0.001);
        assert_abs_diff_eq!(result2.Leq[&desc], 93.98, epsilon = 0.002);
        let lt = result2.Lt.as_ref().expect("Lt should be present");
        let band_lt = lt.get(&desc).expect("descriptor should be in Lt");
        assert!(!band_lt.is_empty(), "Should have time points");
        let time = result2.t.as_ref().expect("time axis should be present");
        assert_eq!(time.len(), band_lt.len());

        Ok(())
    }

    /// Like `test_slm_on_measurement`, but waits for the pre-capture
    /// buffer to fill before starting the recording. The measurement
    /// then contains pre-capture data, which the SLM uses for warm-up.
    #[cfg(feature = "loopback-api")]
    #[test]
    fn test_slm_on_measurement_with_precapture() -> anyhow::Result<()> {
        use crate::config::PRECAP_DURATION;

        let tmpdir = tempfile::tempdir()?;
        let mut smgr = StreamMgr::new_with_devices();

        let in_cfg = DaqConfig {
            api: crate::daq::DaqApiDescriptor::Loopback,
            device_name: "Loopback".into(),
            inchannel_config: vec![
                DaqChannel::defaultAudio("ch0"),
                DaqChannel::defaultAudio("ch1"),
            ],
            outchannel_config: vec![],
            dtype: crate::daq::DataType::F32,
            sampleRateIndex: 0,
            framesPerBlockIndex: 0,
            monitorOutput: false,
        };
        let out_cfg = {
            let mut cfg = in_cfg.clone();
            cfg.inchannel_config = vec![];
            cfg.outchannel_config = vec![
                DaqChannel::defaultAudio("ch0"),
                DaqChannel::defaultAudio("ch1"),
            ];
            cfg
        };

        smgr.startStream(StreamType::Input, &in_cfg)?;
        smgr.startStream(StreamType::Output, &out_cfg)?;
        smgr.setSiggenSource(SourceDescriptor::Sine {
            frequency: (1000.).try_into().unwrap(),
        })?;
        smgr.siggenCommand(SiggenCommand::SetMuteAllChannels { mute: false })?;
        // Set amplitude to √2 → RMS = 1 → 94 dB SPL
        smgr.siggenCommand(SiggenCommand::SetAllGains {
            g: sq2.try_into().unwrap(),
        })?;

        // Wait for the pre-capture buffer to fill (at least 3 seconds).
        std::thread::sleep(PRECAP_DURATION + Duration::from_secs(1));

        let duration = Duration::from_secs(2);
        let rec_settings = RecordSettings::new(
            "slm_precap_test",
            false,
            Some(duration),
            None,
            None,
            Some(MeasurementType::SoundLevels {}),
            false,
            Some(tmpdir.path()),
            vec![],
        )?;

        let mut recording = Recording::new(rec_settings, &mut smgr)?;
        loop {
            if matches!(recording.status(), RecordStatus::Finished { .. }) {
                break;
            }
            std::thread::sleep(Duration::from_millis(10));
        }
        recording.stop();

        let m = recording.getMeasurement().unwrap();
        drop(recording);
        drop(smgr);

        // Verify pre-capture data exists
        let precap_n = m.read().precapture_nframes();
        assert!(precap_n > 0, "Measurement should have pre-capture data");

        let fs = *m.read().samplerate();
        let fsp = StrictlyPositive::new(fs)?;
        let desc = StandardFilterDescriptor::Overall().unwrap();
        let slm_settings = SLMSettingsBuilder::default()
            .fs(fsp)
            .timeWeighting(TimeWeighting::Fast {})
            .freqWeighting(FreqWeighting::Z)
            .filterDescriptors([desc])
            .build()
            .unwrap();

        // With pre-capture warm-up, Lt should be present from t = 0
        let result = m.SLM(slm_settings, 0, None, None, false).unwrap();
        // Measured: Lmax=93.9821, Lpk=96.9897, Leq=93.9796
        assert_abs_diff_eq!(result.Lmax[&desc], 94.0, epsilon = 0.03);
        assert_abs_diff_eq!(result.Lpk[&desc], 96.99, epsilon = 0.001);
        assert_abs_diff_eq!(result.Leq[&desc], 93.98, epsilon = 0.002);
        let lt = result.Lt.as_ref().expect("Lt should be present");
        let band_lt = lt.get(&desc).expect("descriptor should be in Lt");
        assert!(!band_lt.is_empty(), "Should have time points");
        let time = result.t.as_ref().expect("time axis should be present");
        assert_eq!(time.len(), band_lt.len());
        // The first time value should be near zero (pre-capture handled warm-up)
        assert!(time[0] < 0.1, "First output should start near t=0");
        // Measured: last_lt=93.9772
        assert_abs_diff_eq!(*band_lt.last().unwrap(), 93.98, epsilon = 0.005);

        Ok(())
    }
}