measure.hpp

Go to the documentation of this file.

// sbl/dsp/analysis/measure.hpp — Measuring a rendered signal (Cross-cutting — host only)
//
// The measurements that made the bowed string trustworthy, in one place so a
// new model gets them for the cost of a call rather than a copy (AP-041
// Phase 3; RPT-035 Part 4 found the period estimator five times). Host
// tooling and tests only: double arithmetic, std::vector, no ISR promises.
// Never included from audio code.
//
//   bin_energy / bin_power    one frequency's share of a window (Goertzel;
//                             Hann-weighted energy, or rectangular mean power)
//   autocorr_period           the period near an expected lag, interpolated
//   harmonic_profile          level, centroid, each harmonic in dB below the strongest
//   measure_window            periodicity, slips per period, fundamental share,
//                             high-band share, peak, DC, stick fraction
//   magnitude_at              a processor's magnitude at one frequency, by sine drive
//   reflectance               a termination's magnitude and phase delay at one frequency
//   decay_fit                 one partial's decay in dB/s from windowed energies
 
#ifndef SBL_DSP_ANALYSIS_MEASURE_HPP_
#define SBL_DSP_ANALYSIS_MEASURE_HPP_
 
#include <cmath>
#include <cstdint>
#include <vector>
 
namespace sbl::dsp::analysis {
 
inline constexpr double PI = 3.14159265358979323846;
inline constexpr int MAX_HARMONICS = 16;
 
// ─── One frequency's share of a window ───────────────────────────────
 
enum class Window { Rectangular, Hann };
 
/// |Σ w[i] x[i] e^(−jωi)|² by Goertzel's recursion, w the window. The one
/// projection both energy measures are built on.
template<typename T>
inline double projection_power(const T* x, int n, double hz, double fs, Window window) {
    const double w = 2.0 * PI * hz / fs;
    const double coeff = 2.0 * std::cos(w);
    double s1 = 0.0, s2 = 0.0;
    for (int i = 0; i < n; ++i) {
        const double weight = window == Window::Hann ? 0.5 - 0.5 * std::cos(2.0 * PI * i / n) : 1.0;
        const double s0 = weight * x[i] + coeff * s1 - s2;
        s2 = s1;
        s1 = s0;
    }
    return s1 * s1 + s2 * s2 - coeff * s1 * s2;
}
 
/// Hann-weighted energy at hz, normalised by n². The window keeps harmonics
/// apart when bins sit on a measured, not nominal, pitch.
template<typename T>
inline double bin_energy(const T* x, int n, double hz, double fs) {
    return projection_power(x, n, hz, fs, Window::Hann) / (static_cast<double>(n) * n);
}
 
/// Rectangular power at hz as mean power: a unit-amplitude sine at hz reads 0.5,
/// and the bins of a spectrum sum to mean_power(). The DSP suite's measure.
template<typename T>
inline double bin_power(const T* x, int n, double hz, double fs) {
    return 2.0 * projection_power(x, n, hz, fs, Window::Rectangular) / (static_cast<double>(n) * n);
}
 
template<typename T>
inline double mean_power(const T* x, int n) {
    double e = 0.0;
    for (int i = 0; i < n; ++i) e += static_cast<double>(x[i]) * x[i];
    return n > 0 ? e / n : 0.0;
}
 
// ─── Period ──────────────────────────────────────────────────────────
 
struct Period {
    double samples;       ///< the interpolated period
    double periodicity;   ///< autocorrelation at the period over the energy: 1 is periodic there
    int lag;              ///< the integer lag the peak sat on
};
 
/// Autocorrelation near an expected lag: the best integer lag within
/// ±tolerance of it, refined by a parabola through its neighbours. The mean
/// is removed on a copy; every lag is summed over the same span so the peak
/// is not biased toward short lags.
template<typename T>
inline Period autocorr_period(const T* x, int n, double expected, double tolerance) {
    std::vector<double> y(x, x + n);
    double mean = 0.0;
    for (double v : y) mean += v;
    mean /= n > 0 ? n : 1;
    for (double& v : y) v -= mean;
 
    int lo = static_cast<int>(expected * (1.0 - tolerance));
    int hi = static_cast<int>(expected * (1.0 + tolerance)) + 2;
    if (lo < 1) lo = 1;
    if (hi > n - 2) hi = n - 2;
    const int span = n - hi - 1;
    auto corr = [&](int lag) {
        double acc = 0.0;
        for (int i = 0; i < span; ++i) acc += y[i] * y[i + lag];
        return acc;
    };
    Period p{expected, 0.0, lo};
    if (span <= 0 || hi < lo) return p;
 
    double best = corr(lo);
    for (int lag = lo + 1; lag <= hi; ++lag) {
        const double v = corr(lag);
        if (v > best) { best = v; p.lag = lag; }
    }
    const double r0 = corr(0);
    p.periodicity = r0 > 0.0 ? best / r0 : 0.0;
    const double ym1 = corr(p.lag - 1), y0 = best, yp1 = corr(p.lag + 1);
    const double den = ym1 - 2.0 * y0 + yp1;
    p.samples = p.lag + (den != 0.0 ? 0.5 * (ym1 - yp1) / den : 0.0);
    return p;
}
 
// ─── Harmonic profile ────────────────────────────────────────────────
 
struct Profile {
    double rms_db;
    double centroid_hz;
    double db[MAX_HARMONICS];   ///< each harmonic in dB below the strongest; −120 when absent
    int harmonics;
};
 
/// Level, spectral centroid over the harmonics, and each harmonic's level
/// relative to the strongest, on bins at multiples of f0 (Hann-weighted).
template<typename T>
inline Profile harmonic_profile(const T* x, int n, double f0, double fs, int harmonics) {
    Profile pr{};
    if (harmonics > MAX_HARMONICS) harmonics = MAX_HARMONICS;
    pr.harmonics = harmonics;
    const double energy = mean_power(x, n);
    pr.rms_db = energy > 0.0 ? 10.0 * std::log10(energy) : -120.0;
    double amp[MAX_HARMONICS] = {};
    double strongest = 0.0, weighted = 0.0, total = 0.0;
    for (int h = 1; h <= harmonics; ++h) {
        if (h * f0 >= fs * 0.45) break;
        const double e = bin_energy(x, n, h * f0, fs);
        amp[h - 1] = e;
        strongest = strongest > e ? strongest : e;
        weighted += h * f0 * e;
        total += e;
    }
    pr.centroid_hz = total > 0.0 ? weighted / total : 0.0;
    for (int h = 0; h < MAX_HARMONICS; ++h)
        pr.db[h] = (amp[h] > 0.0 && strongest > 0.0) ? 10.0 * std::log10(amp[h] / strongest) : -120.0;
    return pr;
}
 
// ─── A bowed window ──────────────────────────────────────────────────
 
struct Measure {
    double periodicity;       ///< normalised autocorrelation at the played period (1 = periodic there)
    double measured_period;   ///< samples, from the autocorrelation peak near the nominal
    double slips_per_period;  ///< releases per period: 1 is Helmholtz motion, 2 an octave regime
    double h1_percent;        ///< the fundamental's share of the first `harmonics` harmonics' energy
    double hf_percent;        ///< second-difference energy share: the harsh high band
    double peak;
    double dc;
    int stick_percent;
};
 
/// Judge a window of a bowed model's output. The mean is removed from `x`
/// in place, so a profile taken afterwards sees the same signal. Stick and
/// slip counts come from the caller, who watched the exciter while it played.
template<typename T>
inline Measure measure_window(T* x, int n, double nominal_period, double fs, int stick, int slips,
                              int harmonics) {
    Measure m{};
    double mean = 0.0;
    for (int i = 0; i < n; ++i) {
        mean += x[i];
        const double a = std::fabs(static_cast<double>(x[i]));
        m.peak = a > m.peak ? a : m.peak;
    }
    mean /= n > 0 ? n : 1;
    m.dc = mean;
    for (int i = 0; i < n; ++i) x[i] -= static_cast<T>(mean);
    m.stick_percent = n > 0 ? stick * 100 / n : 0;
 
    const Period p = autocorr_period(x, n, nominal_period, 0.03);
    m.periodicity = p.periodicity;
    m.measured_period = p.samples;
    m.slips_per_period = n > 0 ? slips * m.measured_period / n : 0.0;
 
    const double f0 = fs / m.measured_period;
    double total = 0.0, fundamental = 0.0;
    for (int h = 1; h <= harmonics && h * f0 < fs * 0.45; ++h) {
        const double e = bin_energy(x, n, h * f0, fs);
        if (h == 1) fundamental = e;
        total += e;
    }
    m.h1_percent = total > 0.0 ? 100.0 * fundamental / total : 0.0;
 
    const double r0 = mean_power(x, n) * n;
    double hf = 0.0;
    for (int i = 2; i < n; ++i) {
        const double d2 = static_cast<double>(x[i]) - 2.0 * x[i - 1] + x[i - 2];
        hf += d2 * d2 * 0.0625;
    }
    m.hf_percent = r0 > 0.0 ? 100.0 * hf / r0 : 0.0;
    return m;
}
 
/// One release per period, give or take a window edge.
inline bool helmholtz(const Measure& m, double min_slips = 0.9, double max_slips = 1.1) {
    return m.slips_per_period >= min_slips && m.slips_per_period <= max_slips;
}
 
// ─── Probes ──────────────────────────────────────────────────────────
 
/// A processor's magnitude at hz: drive it with a unit sine for `settle`
/// samples, then project its output on the drive over `n`.
template<typename F>
inline double magnitude_at(double hz, double fs, F&& process, int settle = 4000, int n = 20000) {
    const double w = 2.0 * PI * hz / fs;
    double re = 0.0, im = 0.0;
    for (int i = 0; i < settle; ++i) process(static_cast<float>(std::sin(w * i)));
    for (int i = settle; i < settle + n; ++i) {
        const double y = process(static_cast<float>(std::sin(w * i)));
        re += y * std::sin(w * i);
        im += y * std::cos(w * i);
    }
    return std::sqrt(re * re + im * im) * 2.0 / n;
}
 
struct Reflectance {
    double magnitude;
    double phase_delay_samples;   ///< beyond a wall's π: positive when a loop runs longer for it
};
 
/// A termination's reflection for a unit sine arriving at hz.
template<typename F>
inline Reflectance reflectance(double hz, double fs, F&& reflect, int settle = -1, int n = -1) {
    if (settle < 0) settle = static_cast<int>(fs * 0.2);
    if (n < 0) n = static_cast<int>(fs * 0.5);
    const double w = 2.0 * PI * hz / fs;
    for (int i = 0; i < settle; ++i) reflect(static_cast<float>(std::sin(w * i)));
    double re = 0.0, im = 0.0;
    for (int i = settle; i < settle + n; ++i) {
        const double r = reflect(static_cast<float>(std::sin(w * i)));
        re += r * std::cos(w * i);
        im += r * std::sin(w * i);
    }
    const double a = 2.0 * im / n, b = 2.0 * re / n;   // r ≈ a sin + b cos
    Reflectance out{std::sqrt(a * a + b * b), 0.0};
    double extra = std::atan2(b, a) - PI;              // the phase beyond a wall's
    while (extra > PI) extra -= 2.0 * PI;
    while (extra < -PI) extra += 2.0 * PI;
    out.phase_delay_samples = -extra / w;
    return out;
}
 
struct Decay {
    double db_per_s;
    int windows_fit;   ///< windows above the floor that the line was fitted through
};
 
/// One partial's decay: energy at hz in consecutive windows, a straight line
/// through dB against time until 60 dB down or the floor.
template<typename T>
inline Decay decay_fit(const T* x, int n, double hz, int window, double fs, double floor = 1e-12) {
    Decay d{0.0, 0};
    const int windows = window > 0 ? n / window : 0;
    double first_db = 0.0, sx = 0, sy = 0, sxx = 0, sxy = 0;
    for (int w = 0; w < windows; ++w) {
        const double e = bin_energy(x + w * window, window, hz, fs);
        if (e < floor) break;
        const double db = 10.0 * std::log10(e);
        if (w == 0) first_db = db;
        const double t = (w + 0.5) * window / fs;
        sx += t; sy += db; sxx += t * t; sxy += t * db;
        ++d.windows_fit;
        if (d.windows_fit >= 2 && db < first_db - 60.0) break;
    }
    const int k = d.windows_fit;
    d.db_per_s = k >= 2 ? (k * sxy - sx * sy) / (k * sxx - sx * sx) : 0.0;
    return d;
}
 
} // namespace sbl::dsp::analysis
 
#endif // SBL_DSP_ANALYSIS_MEASURE_HPP_