From 57d4c5569fc1014458133c730aa399b9db1ae3b4 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Wed, 29 Jul 2026 15:14:28 +0000 Subject: [PATCH 1/3] add metrics for libraries --- README.md | 4 +- ROADMAP.md | 19 -- doc/analysis/Decibels.md | 98 +++++++++++ doc/analysis/README.md | 1 + doc/math/README.md | 5 +- doc/math/StepResponseMetrics.md | 107 ++++++++++++ numerical/analysis/CMakeLists.txt | 2 + numerical/analysis/Decibels.cpp | 9 + numerical/analysis/Decibels.hpp | 57 ++++++ numerical/analysis/test/CMakeLists.txt | 1 + numerical/analysis/test/TestDecibels.cpp | 60 +++++++ numerical/math/CMakeLists.txt | 2 + numerical/math/StepResponseMetrics.cpp | 10 ++ numerical/math/StepResponseMetrics.hpp | 124 +++++++++++++ numerical/math/test/CMakeLists.txt | 1 + .../math/test/TestStepResponseMetrics.cpp | 164 ++++++++++++++++++ 16 files changed, 641 insertions(+), 23 deletions(-) create mode 100644 doc/analysis/Decibels.md create mode 100644 doc/math/StepResponseMetrics.md create mode 100644 numerical/analysis/Decibels.cpp create mode 100644 numerical/analysis/Decibels.hpp create mode 100644 numerical/analysis/test/TestDecibels.cpp create mode 100644 numerical/math/StepResponseMetrics.cpp create mode 100644 numerical/math/StepResponseMetrics.hpp create mode 100644 numerical/math/test/TestStepResponseMetrics.cpp diff --git a/README.md b/README.md index 28eb51c2..db214f76 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal | Category | Description | |--------------------------------------------------------------------|----------------------------------------------------------------------| -| [Analysis](doc/analysis/README.md) | FFT, Real-Input FFT (RFFT), Power Spectral Density, DCT, Window Functions, Signal Detectors, Convolution & Correlation, Goertzel Algorithm | +| [Analysis](doc/analysis/README.md) | FFT, Real-Input FFT (RFFT), Power Spectral Density, DCT, Window Functions, Signal Detectors, Convolution & Correlation, Goertzel Algorithm, Decibels | | [Control Analysis](doc/control_analysis/README.md) | Frequency Response, Root Locus, Controllability/Observability Matrices & Gramians | | [Controllers](doc/controllers/README.md) | Bang-Bang/Hysteresis, PID, LQR, LQI (Integral/Servo State Feedback), MPC, Saturation, Rate Limiter, Slew-Limited Saturation, Feedforward/2-DOF, Gain-Scheduled Controller, Lead-Lag Compensator, Luenberger Observer | | [Dynamics](doc/dynamics/README.md) | Euler-Lagrange, Newton-Euler, Recursive Newton-Euler, ABA | @@ -26,7 +26,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal | [Neural Network](doc/neural_network/README.md) | Layers, activations, losses, model | | [Optimization](doc/optimization/README.md) | Gradient Descent | | [Regularization](doc/regularization/README.md) | L1 (Lasso), L2 (Ridge) | -| [Math](doc/math/README.md) | CORDIC, Quaternion | +| [Math](doc/math/README.md) | CORDIC, Quaternion, Step Response Metrics | | [Solvers](doc/solvers/README.md) | Gaussian Elimination, Levinson-Durbin, Durand-Kerner, Cholesky, DARE, Runge-Kutta ODE Integrators (RK4 + Dormand-Prince) | | [Performance Optimization](doc/performance-optimization/README.md) | Compiler optimizations, SIMD | diff --git a/ROADMAP.md b/ROADMAP.md index 554ed89b..a4d39ed8 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -48,8 +48,6 @@ Difficulty legend: | 45 | IIR filter design (Butterworth/Chebyshev + bilinear) | `filters/passive` | ★★★★★ | | 46 | H∞ state-feedback control | `robust_control` (new) | ★★★★★ | | 47 | Model Reference Adaptive Control (MRAC) | `nonlinear_control` (new) | ★★★★★ | -| 48 | Decibel & magnitude-response helpers | `analysis` | ★☆☆☆☆ | -| 49 | Step / transient-response metrics | `math` | ★★☆☆☆ | | 50 | Matrix norms & condition number | `math` | ★★★☆☆ | | 51 | Spectral radius / discrete stability margin | `math` | ★★★☆☆ | | 52 | Estimator consistency metrics (NEES / NIS) | `estimators` | ★★★☆☆ | @@ -561,23 +559,6 @@ and [`control_analysis/FrequencyResponse`](numerical/control_analysis/FrequencyR magnitude/phase. The items below are the missing pieces. All are **float-only**, no-heap, and operate on bounded `math::Vector`/`math::Matrix` inputs; tests are `TEST_F` on `float`. -### 48. Decibel & magnitude-response helpers ★☆☆☆☆ — `analysis` -- **What:** `ToDecibels(ratio)` / `FromDecibels`, plus convenience magnitude(-in-dB) and attenuation - helpers layered over `control_analysis::FrequencyResponse`. -- **Metric value:** M2 (frequency response) — pass-band ripple / stop-band attenuation in dB, the - natural unit for filter and controller tests. -- **Algorithm:** `20·log10(·)`; guard the zero/`-inf` case with a floor. -- **Reuses:** `math::TrigonometricFunctions`/`std::log10`, existing `FrequencyResponse`. - -### 49. Step / transient-response metrics ★★☆☆☆ — `math` -- **What:** From a bounded step-response `Vector`: `RiseTime` (10–90 %), `SettlingTime` (±band), - `PercentOvershoot`, `PeakTime`, `SteadyStateError`. -- **Metric value:** M3 (time-response) — the core acceptance criteria for every controller - (`controllers/`), `LinearTimeInvariant`, and IIR filter. -- **Algorithm:** single forward pass over the sampled response against the reference/steady value; - standard control-systems definitions. -- **Reuses:** `math::Vector`, `math::Statistics` for the steady-state estimate. - ### 50. Matrix norms & condition number ★★★☆☆ — `math` - **What:** `FrobeniusNorm`, `OneNorm`, `InfinityNorm` on `Matrix`; `Vector` `Norm`/`Normalize`; `ConditionNumber` estimate. diff --git a/doc/analysis/Decibels.md b/doc/analysis/Decibels.md new file mode 100644 index 00000000..08b23722 --- /dev/null +++ b/doc/analysis/Decibels.md @@ -0,0 +1,98 @@ +# Decibels & Magnitude-Response Helpers + +## Overview & Motivation + +Audio, filter, and control-system specifications express signal levels and filter performance in decibels (dB) because the human auditory system and most engineering metrics scale logarithmically with amplitude ratio. Specifying a stop-band attenuation of 60 dB or a pass-band ripple of 0.1 dB is natural and compact; the equivalent linear ratios (1 000 : 1 and 1.01161 : 1) are not. A small set of conversion primitives — `ToDecibels`, `FromDecibels`, and two derived helpers for attenuation and ripple — centralises this conversion and eliminates scattered, error-prone inline `20·log10` expressions throughout the rest of the library. + +## Mathematical Theory + +### Magnitude Decibel Conversion + +For a positive amplitude ratio $r > 0$, the equivalent level in decibels is: + +$$L_{\mathrm{dB}} = 20 \log_{10}(r)$$ + +The factor 20 (rather than 10) reflects the voltage/pressure convention: power is proportional to the square of amplitude, so a doubling of amplitude ($r = 2$) gives a 6 dB increase, matching the $10 \log_{10}(4) = 6.02$ dB power equivalent. + +### Inverse Conversion + +$$r = 10^{L_{\mathrm{dB}}/20}$$ + +This inverse is exact for all finite $L_{\mathrm{dB}}$; no guard is needed on the output side. + +### Zero and Negative Input Guard + +$\log_{10}(0) = -\infty$; negative ratios are physically meaningless. Both cases are mapped to a finite floor value $L_{\min}$ chosen well below any engineering specification of interest: + +$$L_{\mathrm{dB}} = \max\!\left(20\log_{10}(r),\; L_{\min}\right), \quad r > 0$$ +$$L_{\mathrm{dB}} = L_{\min}, \quad r \leq 0$$ + +A floor of $-160\,\mathrm{dB}$ corresponds to an amplitude ratio below $10^{-8}$, safely beyond the dynamic range of any practical floating-point computation in 32-bit single precision. + +### Derived Helpers + +**Stop-band attenuation** between a pass-band ratio $r_p$ and a stop-band ratio $r_s$: + +$$A = L_{\mathrm{dB}}(r_p) - L_{\mathrm{dB}}(r_s)$$ + +**Pass-band ripple** between the maximum and minimum in-band ratios $r_{\max}$ and $r_{\min}$: + +$$\Delta = L_{\mathrm{dB}}(r_{\max}) - L_{\mathrm{dB}}(r_{\min})$$ + +Both are simple differences in decibel space, exploiting the logarithm identity $\log(a/b) = \log a - \log b$. + +## Complexity Analysis + +| Operation | Time | Space | Notes | +|-----------------|----------|-------|------------------------------------| +| `ToDecibels` | O(1) | O(1) | One `log10` + one `max` + one `mul` | +| `FromDecibels` | O(1) | O(1) | One `pow` | +| `AttenuationDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction | +| `RippleDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction | + +No state, no buffers. All operations are pure functions. + +## Step-by-Step Walkthrough + +Converting a ratio of 10 to decibels: + +1. Input $r = 10$; guard passes ($r > 0$). +2. Compute $20 \cdot \log_{10}(10) = 20 \cdot 1 = 20$. +3. Apply floor: $\max(20, -160) = 20$. +4. Output: $20\,\mathrm{dB}$. + +Round-trip for $r = 0.5$: + +1. `ToDecibels(0.5)` = $20 \cdot \log_{10}(0.5) \approx -6.0206\,\mathrm{dB}$. +2. `FromDecibels(-6.0206)` = $10^{-6.0206/20} \approx 0.5$. + +## Pitfalls & Edge Cases + +- Passing $r = 0$ produces $-\infty$ from `log10`; the floor guard prevents propagation into downstream computations. +- Negative ratios indicate a programming error (signed sample values must not be passed directly as ratios without taking absolute value first); they are silently floored rather than raising an exception, consistent with the no-exception policy. +- With `fast-math` enabled, the compiler may fuse or reorder floating-point operations. The `log10` result is still monotone and the floor remains correct because it uses `std::max`, which is not reordered away. +- `FromDecibels` has no floor; at very large positive dB values the result overflows to `+inf` in float — this is expected behaviour for out-of-range inputs. + +## Variants & Generalizations + +- Power decibels use $10 \log_{10}(\cdot)$ (factor 10 rather than 20). The amplitude convention used here ($\times 20$) is correct for voltage, pressure, and filter transfer-function magnitude. +- Field-quantity vs. power-quantity disambiguation: IEEE 60268 / IEC 61672 mandate $20 \log_{10}$ for sound pressure level; the same convention applies to filter magnitude response. +- The floor can be parameterised if a stricter or looser sentinel is required; the default of $-160\,\mathrm{dB}$ is conservative for 32-bit float. + +## Applications + +- Filter specification: stop-band attenuation and pass-band ripple in dB are the primary acceptance criteria for IIR/FIR designs. +- Frequency response plots: `FrequencyResponse::Calculate()` already returns magnitude in dB using $20 \log_{10}$; these helpers provide the same conversion for ad-hoc analysis. +- Controller gain margin is expressed in dB; converting from a linear ratio with `ToDecibels` avoids duplication. +- Audio dynamic processing (compressor thresholds, limiter ceilings) and acoustic measurement both use dB natively. + +## Connections to Other Algorithms + +- `control_analysis::FrequencyResponse` internally applies $20 \log_{10}(\|H\|)$ on its magnitude output vector; these helpers are the scalar equivalent exposed for library consumers. +- Pass-band ripple computed by `RippleDb` feeds directly into filter-design acceptance testing alongside the step/transient-response metrics in the evaluation primitives family. + +## References & Further Reading + +- Proakis, J. & Manolakis, D., "Digital Signal Processing", 4th ed., Prentice Hall, 2007 — Appendix A (decibel notation). +- Zolzer, U., "DAFX: Digital Audio Effects", 2nd ed., Wiley, 2011 — Chapter 2 (level and gain in dB). +- IEC 61672-1:2013, "Electroacoustics — Sound level meters — Part 1: Specifications." diff --git a/doc/analysis/README.md b/doc/analysis/README.md index 596d7c31..f55564aa 100644 --- a/doc/analysis/README.md +++ b/doc/analysis/README.md @@ -11,6 +11,7 @@ Signal analysis algorithms for frequency-domain decomposition and spectral estim | [Power Spectral Density](PowerDensitySpectrum.md) | Estimation of signal power distribution across frequencies using Welch's method | | [Discrete Cosine Transform](DiscreteCosineTransform.md) | Real-valued frequency decomposition via cosine basis functions, computed through FFT | | [Signal Detectors](SignalDetectors.md) | Peak hold, zero-crossing counter, and RMS envelope detectors for real-time signal monitoring | +| [Decibels](Decibels.md) | `ToDecibels` / `FromDecibels` conversion helpers with zero-floor guard, plus attenuation and ripple utilities | | [Goertzel Algorithm](GoertzelAlgorithm.md) | Single-bin DFT via a second-order recurrence for O(N) tone detection with O(1) memory | ## Sub-domains diff --git a/doc/math/README.md b/doc/math/README.md index 4854645f..e5e1c108 100644 --- a/doc/math/README.md +++ b/doc/math/README.md @@ -6,5 +6,6 @@ Core mathematical primitives for numerical computation. | Algorithm | Description | |-----------------------------|-----------------------------------------------------------------------------------------------| -| [CORDIC](Cordic.md) | Iterative shift-add engine for sin/cos, atan2, magnitude, and vector rotation — no multiplier | -| [Quaternion](Quaternion.md) | Unit-quaternion rotation type: Hamilton product, SLERP, rotation-matrix and Euler conversions | +| [CORDIC](Cordic.md) | Iterative shift-add engine for sin/cos, atan2, magnitude, and vector rotation — no multiplier | +| [Quaternion](Quaternion.md) | Unit-quaternion rotation type: Hamilton product, SLERP, rotation-matrix and Euler conversions | +| [Step Response Metrics](StepResponseMetrics.md) | Rise time, settling time, percent overshoot, peak time, and steady-state error from a bounded step-response vector | diff --git a/doc/math/StepResponseMetrics.md b/doc/math/StepResponseMetrics.md new file mode 100644 index 00000000..273ea749 --- /dev/null +++ b/doc/math/StepResponseMetrics.md @@ -0,0 +1,107 @@ +# Step / Transient-Response Metrics + +## Overview & Motivation + +When a control system or filter receives a step input, its output traces a transient trajectory before settling at the final value. Quantifying that trajectory with standardised scalar metrics — rise time, settling time, percent overshoot, peak time, and steady-state error — is the primary acceptance test for any closed-loop design. These metrics translate the raw sample sequence into the language of control specifications, allowing automated pass/fail decisions without manual inspection of time-domain plots. + +## Mathematical Theory + +### Definitions + +Let $y[k]$, $k = 0, \ldots, N-1$ be the sampled step response and $y_{ss}$ the steady-state value. The sample period is $\Delta t$. + +**Rise Time** $T_r$ + +The elapsed time for the response to travel from 10 % to 90 % of steady state: + +$$T_r = (k_{90} - k_{10})\,\Delta t$$ + +where $k_{10} = \min\{k : y[k] \ge 0.1\,y_{ss}\}$ and $k_{90} = \min\{k \ge k_{10} : y[k] \ge 0.9\,y_{ss}\}$. + +**Settling Time** $T_s$ + +The first time after which the response remains permanently inside the band $[(1-\delta)y_{ss},\,(1+\delta)y_{ss}]$ (typically $\delta = 0.02$): + +$$T_s = (k^* + 1)\,\Delta t, \quad k^* = \max\{k : |y[k] - y_{ss}| > \delta\,|y_{ss}|\}$$ + +**Percent Overshoot** $\%OS$ + +$$\%OS = 100\,\frac{y_{\max} - y_{ss}}{y_{ss}}, \quad y_{\max} = \max_k y[k]$$ + +For an underdamped second-order system with damping ratio $\zeta$: + +$$\%OS = 100\,\exp\!\left(-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}\right)$$ + +**Peak Time** $T_p$ + +$$T_p = k_p\,\Delta t, \quad k_p = \arg\max_k y[k]$$ + +For a continuous underdamped second-order system with natural frequency $\omega_n$: + +$$T_p = \frac{\pi}{\omega_n\sqrt{1-\zeta^2}}$$ + +**Steady-State Error** $e_{ss}$ + +$$e_{ss} = r - \bar{y}_{\text{tail}}$$ + +where $r$ is the reference (command) value and $\bar{y}_{\text{tail}}$ is the mean of the final quarter of the response buffer, providing a robust estimate of the achieved steady state. + +## Complexity Analysis + +| Case | Time | Space | Notes | +|---------|----------|--------|--------------------------------------------| +| All | $O(N)$ | $O(1)$ | Single forward pass; no auxiliary storage | + +Each metric requires at most one traversal of the $N$-element vector. The tail-mean for steady-state error adds a constant-fraction second scan of the same data — still $O(N)$ total. + +## Step-by-Step Walkthrough + +Consider a 10-sample ramp to $y_{ss} = 1$ followed by a constant plateau (N = 20): + +``` +k: 0 1 2 3 4 5 6 7 8 9 10 11 … +y: 0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1 1 … +``` + +- **Rise Time:** $k_{10} = 1$ (first sample $\ge 0.1$), $k_{90} = 9$ (first sample $\ge 0.9$). $T_r = 8\,\Delta t$. +- **Settling Time:** With $\delta = 0.02$, last sample outside the band is $k = 9$. $T_s = 10\,\Delta t$. +- **Percent Overshoot:** $y_{\max} = 1.0 = y_{ss}$, so $\%OS = 0$. +- **Peak Time:** $k_p = 10$ (first occurrence of max). $T_p = 10\,\Delta t$. +- **Steady-State Error:** Tail mean $= 1.0$, reference $= 1.0$. $e_{ss} = 0$. + +## Pitfalls & Edge Cases + +**Zero steady state.** Division by $y_{ss}$ in percent overshoot is guarded; the function returns zero when $y_{ss} = 0$ to avoid a NaN. + +**Non-monotone ramp.** If the response crosses 90 % before 10 % (e.g., DC offset or wrong initial condition), $k_{10}$ may be found after the first 90 % crossing. The implementation returns the first pair that satisfies the threshold order. + +**Oscillatory settling.** Settling time is defined as the last time the trajectory leaves the band, not the first time it enters it. Repeated crossings near the boundary extend the metric correctly. + +**Finite buffer.** With a bounded vector of length $N$, if the response has not yet settled by the final sample, `SettlingTime` returns $N\,\Delta t$ and `RiseTime` returns $(N-1)\,\Delta t$ as conservative bounds. + +**Tail-mean length.** Using the last $\lfloor N/4 \rfloor + 1$ samples for the steady-state estimate assumes the transient has decayed to within numerical noise by that point. Poorly chosen $N$ relative to the system time constant degrades the estimate. + +## Variants & Generalizations + +- **Delay Time** $T_d$: the time to reach 50 % of steady state — obtainable with the same threshold-scan pattern. +- **Band-relative rise time**: using a band other than 10–90 % (e.g., 20–80 %) is a trivial parameter change. +- **Multi-channel:** applying the scalar functions element-wise to each row of a response matrix generalises to MIMO systems without algorithmic change. + +## Applications + +- Automated controller tuning acceptance: verify that a PID or LQR design meets specification ($T_r < T_{r,\text{spec}}$, $\%OS < \%OS_{\text{spec}}$, etc.). +- Filter characterisation: measure the transient of a step fed through an IIR or FIR filter. +- Hardware-in-the-loop test harnesses: compute metrics directly from sampled actuator responses. + +## Connections to Other Algorithms + +- **Statistics** (this library): the tail-mean for steady-state error replicates the `Mean` function on a sub-range. +- **LinearTimeInvariant**: the primary source of step responses whose metrics are evaluated here. +- **Filters/active** (Kalman, EKF): step-excitation tests use these metrics to validate estimator transient behaviour. +- **Controllers**: PID and LQR tuning loops iterate until all five metrics satisfy design targets. + +## References & Further Reading + +- K. J. Åström and R. M. Murray, *Feedback Systems: An Introduction for Scientists and Engineers*, Princeton University Press, 2008. Chapter 10. +- G. F. Franklin, J. D. Powell, and A. Emami-Naeini, *Feedback Control of Dynamic Systems*, 8th ed., Pearson, 2019. Chapter 3. +- N. S. Nise, *Control Systems Engineering*, 8th ed., Wiley, 2019. Chapter 4. diff --git a/numerical/analysis/CMakeLists.txt b/numerical/analysis/CMakeLists.txt index 7b46ff3f..13348ade 100644 --- a/numerical/analysis/CMakeLists.txt +++ b/numerical/analysis/CMakeLists.txt @@ -13,6 +13,7 @@ target_link_libraries(numerical.analysis ${NUMERICAL_VISIBILITY} target_sources(numerical.analysis PRIVATE ConvolutionCorrelation.hpp + Decibels.hpp DiscreteCosineTransform.hpp FastFourierTransform.hpp FastFourierTransformRadix2Impl.hpp @@ -24,6 +25,7 @@ target_sources(numerical.analysis PRIVATE numerical_add_coverage_sources(numerical.analysis ConvolutionCorrelation.cpp + Decibels.cpp FastFourierTransformRadix2Impl.cpp GoertzelAlgorithm.cpp PowerDensitySpectrum.cpp diff --git a/numerical/analysis/Decibels.cpp b/numerical/analysis/Decibels.cpp new file mode 100644 index 00000000..78cff0bd --- /dev/null +++ b/numerical/analysis/Decibels.cpp @@ -0,0 +1,9 @@ +#include "numerical/analysis/Decibels.hpp" + +namespace analysis +{ + template float ToDecibels(float); + template float FromDecibels(float); + template float AttenuationDb(float, float); + template float RippleDb(float, float); +} diff --git a/numerical/analysis/Decibels.hpp b/numerical/analysis/Decibels.hpp new file mode 100644 index 00000000..86cdea98 --- /dev/null +++ b/numerical/analysis/Decibels.hpp @@ -0,0 +1,57 @@ +#pragma once + +#if defined(__GNUC__) || defined(__clang__) +#pragma GCC optimize("O3", "fast-math") +#endif + +#include "numerical/math/CompilerOptimizations.hpp" +#include +#include +#include + +namespace analysis +{ + template + struct DecibelFloor + { + static_assert(std::is_floating_point_v, "DecibelFloor supports floating-point types only"); + static constexpr T value{ T{ -160 } }; + }; + + template + OPTIMIZE_FOR_SPEED T ToDecibels(T ratio) + { + static_assert(std::is_floating_point_v, "ToDecibels supports floating-point types only"); + if (ratio <= T{ 0 }) + return DecibelFloor::value; + return std::max(T{ 20 } * std::log10(ratio), DecibelFloor::value); + } + + template + OPTIMIZE_FOR_SPEED T FromDecibels(T db) + { + static_assert(std::is_floating_point_v, "FromDecibels supports floating-point types only"); + return std::pow(T{ 10 }, db / T{ 20 }); + } + + template + OPTIMIZE_FOR_SPEED T AttenuationDb(T passbandRatio, T stopbandRatio) + { + static_assert(std::is_floating_point_v, "AttenuationDb supports floating-point types only"); + return ToDecibels(passbandRatio) - ToDecibels(stopbandRatio); + } + + template + OPTIMIZE_FOR_SPEED T RippleDb(T maxRatio, T minRatio) + { + static_assert(std::is_floating_point_v, "RippleDb supports floating-point types only"); + return ToDecibels(maxRatio) - ToDecibels(minRatio); + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template float ToDecibels(float); + extern template float FromDecibels(float); + extern template float AttenuationDb(float, float); + extern template float RippleDb(float, float); +#endif +} diff --git a/numerical/analysis/test/CMakeLists.txt b/numerical/analysis/test/CMakeLists.txt index de8a4acb..8458ad3b 100644 --- a/numerical/analysis/test/CMakeLists.txt +++ b/numerical/analysis/test/CMakeLists.txt @@ -10,6 +10,7 @@ target_link_libraries(numerical.analysis_test PUBLIC target_sources(numerical.analysis_test PRIVATE TestConvolutionCorrelation.cpp + TestDecibels.cpp TestDiscreteCosineTransform.cpp TestFastFourierTransformRadix2Impl.cpp TestGoertzelAlgorithm.cpp diff --git a/numerical/analysis/test/TestDecibels.cpp b/numerical/analysis/test/TestDecibels.cpp new file mode 100644 index 00000000..4c02ec31 --- /dev/null +++ b/numerical/analysis/test/TestDecibels.cpp @@ -0,0 +1,60 @@ +#include "numerical/analysis/Decibels.hpp" +#include "numerical/math/Tolerance.hpp" +#include "gmock/gmock.h" +#include + +namespace +{ + class TestDecibels : public ::testing::Test + {}; +} + +TEST_F(TestDecibels, ratio_ten_yields_twenty_db) +{ + EXPECT_NEAR(analysis::ToDecibels(10.0f), 20.0f, math::Tolerance()); +} + +TEST_F(TestDecibels, ratio_hundred_yields_forty_db) +{ + EXPECT_NEAR(analysis::ToDecibels(100.0f), 40.0f, math::Tolerance()); +} + +TEST_F(TestDecibels, ratio_half_yields_minus_six_db) +{ + EXPECT_NEAR(analysis::ToDecibels(0.5f), -6.0206f, 1e-3f); +} + +TEST_F(TestDecibels, ratio_one_yields_zero_db) +{ + EXPECT_NEAR(analysis::ToDecibels(1.0f), 0.0f, math::Tolerance()); +} + +TEST_F(TestDecibels, zero_ratio_returns_floor) +{ + EXPECT_NEAR(analysis::ToDecibels(0.0f), analysis::DecibelFloor::value, math::Tolerance()); +} + +TEST_F(TestDecibels, negative_ratio_returns_floor) +{ + EXPECT_NEAR(analysis::ToDecibels(-1.0f), analysis::DecibelFloor::value, math::Tolerance()); +} + +TEST_F(TestDecibels, from_decibels_twenty_returns_ten) +{ + EXPECT_NEAR(analysis::FromDecibels(20.0f), 10.0f, math::Tolerance()); +} + +TEST_F(TestDecibels, round_trip_preserves_ratio) +{ + EXPECT_NEAR(analysis::FromDecibels(analysis::ToDecibels(0.5f)), 0.5f, math::Tolerance()); +} + +TEST_F(TestDecibels, attenuation_db_computes_difference) +{ + EXPECT_NEAR(analysis::AttenuationDb(1.0f, 0.01f), 40.0f, math::Tolerance()); +} + +TEST_F(TestDecibels, ripple_db_computes_passband_variation) +{ + EXPECT_NEAR(analysis::RippleDb(1.0f, 0.9f), analysis::ToDecibels(1.0f) - analysis::ToDecibels(0.9f), math::Tolerance()); +} diff --git a/numerical/math/CMakeLists.txt b/numerical/math/CMakeLists.txt index 7d29bb4e..9b10d554 100644 --- a/numerical/math/CMakeLists.txt +++ b/numerical/math/CMakeLists.txt @@ -22,6 +22,7 @@ target_sources(numerical.math PRIVATE RecursiveBuffer.hpp SingleInstructionMultipleData.hpp Statistics.hpp + StepResponseMetrics.hpp Toeplitz.hpp Tolerance.hpp TrigonometricFunctions.hpp @@ -34,6 +35,7 @@ numerical_add_coverage_sources(numerical.math Matrix.cpp QNumber.cpp Quaternion.cpp + StepResponseMetrics.cpp ) add_subdirectory(test) diff --git a/numerical/math/StepResponseMetrics.cpp b/numerical/math/StepResponseMetrics.cpp new file mode 100644 index 00000000..8ac9f69a --- /dev/null +++ b/numerical/math/StepResponseMetrics.cpp @@ -0,0 +1,10 @@ +#include "numerical/math/StepResponseMetrics.hpp" + +namespace math +{ + template float RiseTime(const Vector&, float, float); + template float SettlingTime(const Vector&, float, float, float); + template float PercentOvershoot(const Vector&, float); + template float PeakTime(const Vector&, float); + template float SteadyStateError(const Vector&, float); +} diff --git a/numerical/math/StepResponseMetrics.hpp b/numerical/math/StepResponseMetrics.hpp new file mode 100644 index 00000000..3cc0b0ec --- /dev/null +++ b/numerical/math/StepResponseMetrics.hpp @@ -0,0 +1,124 @@ +#pragma once + +#if defined(__GNUC__) || defined(__clang__) +#pragma GCC optimize("O3", "fast-math") +#endif + +#include "numerical/math/CompilerOptimizations.hpp" +#include "numerical/math/Matrix.hpp" +#include +#include +#include +#include + +namespace math +{ + template + [[nodiscard]] OPTIMIZE_FOR_SPEED T RiseTime(const Vector& response, T steady, T dt = T{ 1 }) + { + static_assert(std::is_floating_point_v, "RiseTime supports floating-point types only"); + + const T lo{ T{ 0.1 } * steady }; + const T hi{ T{ 0.9 } * steady }; + + std::size_t iLo{ 0 }; + bool foundLo{ false }; + + for (std::size_t i = 0; i < Size; ++i) + { + const T val{ response.at(i, 0) }; + if (!foundLo && val >= lo) + { + iLo = i; + foundLo = true; + } + if (foundLo && val >= hi) + return static_cast(i - iLo) * dt; + } + + return static_cast(Size - 1) * dt; + } + + template + [[nodiscard]] OPTIMIZE_FOR_SPEED T SettlingTime(const Vector& response, T steady, T band = T{ 0.02 }, T dt = T{ 1 }) + { + static_assert(std::is_floating_point_v, "SettlingTime supports floating-point types only"); + + const T absThreshold{ std::abs(steady) * band }; + std::size_t lastOutside{ 0 }; + bool anyOutside{ false }; + + for (std::size_t i = 0; i < Size; ++i) + { + if (std::abs(response.at(i, 0) - steady) > absThreshold) + { + lastOutside = i; + anyOutside = true; + } + } + + if (!anyOutside) + return T{ 0 }; + + return static_cast(lastOutside + 1) * dt; + } + + template + [[nodiscard]] OPTIMIZE_FOR_SPEED T PercentOvershoot(const Vector& response, T steady) + { + static_assert(std::is_floating_point_v, "PercentOvershoot supports floating-point types only"); + + if (steady == T{ 0 }) + return T{ 0 }; + + T peak{ response.at(0, 0) }; + for (std::size_t i = 1; i < Size; ++i) + peak = std::max(peak, response.at(i, 0)); + + return T{ 100 } * (peak - steady) / steady; + } + + template + [[nodiscard]] OPTIMIZE_FOR_SPEED T PeakTime(const Vector& response, T dt = T{ 1 }) + { + static_assert(std::is_floating_point_v, "PeakTime supports floating-point types only"); + + std::size_t peakIdx{ 0 }; + T peakVal{ response.at(0, 0) }; + + for (std::size_t i = 1; i < Size; ++i) + { + const T val{ response.at(i, 0) }; + if (val > peakVal) + { + peakVal = val; + peakIdx = i; + } + } + + return static_cast(peakIdx) * dt; + } + + template + [[nodiscard]] OPTIMIZE_FOR_SPEED T SteadyStateError(const Vector& response, T reference) + { + static_assert(std::is_floating_point_v, "SteadyStateError supports floating-point types only"); + + T sum{ T{ 0 } }; + constexpr std::size_t tailStart{ Size - TailSize }; + + for (std::size_t i = tailStart; i < Size; ++i) + sum += response.at(i, 0); + + const T tailMean{ sum / static_cast(TailSize) }; + return reference - tailMean; + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template float RiseTime(const Vector&, float, float); + extern template float SettlingTime(const Vector&, float, float, float); + extern template float PercentOvershoot(const Vector&, float); + extern template float PeakTime(const Vector&, float); + extern template float SteadyStateError(const Vector&, float); +#endif +} diff --git a/numerical/math/test/CMakeLists.txt b/numerical/math/test/CMakeLists.txt index c6c8cc53..55577c1a 100644 --- a/numerical/math/test/CMakeLists.txt +++ b/numerical/math/test/CMakeLists.txt @@ -16,5 +16,6 @@ target_sources(numerical.math_test PRIVATE TestQuaternion.cpp TestRecursiveBuffer.cpp TestStatistics.cpp + TestStepResponseMetrics.cpp TestToeplitz.cpp ) diff --git a/numerical/math/test/TestStepResponseMetrics.cpp b/numerical/math/test/TestStepResponseMetrics.cpp new file mode 100644 index 00000000..417a79b3 --- /dev/null +++ b/numerical/math/test/TestStepResponseMetrics.cpp @@ -0,0 +1,164 @@ +#include "numerical/math/StepResponseMetrics.hpp" +#include "numerical/math/Tolerance.hpp" +#include "gmock/gmock.h" +#include + +namespace +{ + static constexpr std::size_t N = 64; + using Vec = math::Vector; + + class TestStepResponseMetrics : public ::testing::Test + { + protected: + Vec MakeRampPlateau(float steady, std::size_t rampEnd) const + { + Vec v; + for (std::size_t i = 0; i < N; ++i) + { + if (i <= rampEnd) + v.at(i, 0) = steady * static_cast(i) / static_cast(rampEnd); + else + v.at(i, 0) = steady; + } + return v; + } + + Vec MakeUnderdamped(float steady, float overshootFrac, std::size_t peakIdx) const + { + Vec v; + const float peak{ steady * (1.0f + overshootFrac) }; + for (std::size_t i = 0; i < N; ++i) + { + const float t{ static_cast(i) / static_cast(peakIdx) }; + const float envelope{ 1.0f - std::exp(-2.0f * t) }; + const float osc{ overshootFrac * std::exp(-2.0f * t) * std::cos(std::numbers::pi_v * t) }; // NOLINT + v.at(i, 0) = steady * (envelope + osc); + } + (void)peak; + return v; + } + }; +} + +TEST_F(TestStepResponseMetrics, rise_time_ramp_plateau_in_samples) +{ + const float steady{ 1.0f }; + const std::size_t rampEnd{ 20 }; + Vec v{ MakeRampPlateau(steady, rampEnd) }; + + const float result{ math::RiseTime(v, steady) }; + + EXPECT_NEAR(result, 16.0f, 2.0f); +} + +TEST_F(TestStepResponseMetrics, rise_time_with_dt_scales_to_real_time) +{ + const float steady{ 1.0f }; + Vec v{ MakeRampPlateau(steady, 20) }; + const float dt{ 0.01f }; + + const float sampleResult{ math::RiseTime(v, steady) }; + const float timeResult{ math::RiseTime(v, steady, dt) }; + + EXPECT_NEAR(timeResult, sampleResult * dt, math::Tolerance()); +} + +TEST_F(TestStepResponseMetrics, settling_time_at_steady_state_returns_zero) +{ + Vec v; + for (std::size_t i = 0; i < N; ++i) + v.at(i, 0) = 1.0f; + + const float result{ math::SettlingTime(v, 1.0f) }; + + EXPECT_NEAR(result, 0.0f, math::Tolerance()); +} + +TEST_F(TestStepResponseMetrics, settling_time_ramp_plateau_within_band) +{ + const float steady{ 1.0f }; + Vec v{ MakeRampPlateau(steady, 20) }; + + const float result{ math::SettlingTime(v, steady, 0.02f) }; + + EXPECT_GT(result, 0.0f); + EXPECT_LT(result, static_cast(N)); +} + +TEST_F(TestStepResponseMetrics, percent_overshoot_no_overshoot_returns_non_positive) +{ + Vec v{ MakeRampPlateau(1.0f, 30) }; + + const float result{ math::PercentOvershoot(v, 1.0f) }; + + EXPECT_LE(result, 0.0f); +} + +TEST_F(TestStepResponseMetrics, percent_overshoot_known_peak) +{ + Vec v; + for (std::size_t i = 0; i < N; ++i) + v.at(i, 0) = 1.0f; + v.at(10, 0) = 1.2f; + + const float result{ math::PercentOvershoot(v, 1.0f) }; + + EXPECT_NEAR(result, 20.0f, 1e-3f); +} + +TEST_F(TestStepResponseMetrics, percent_overshoot_zero_steady_returns_zero) +{ + Vec v{ MakeRampPlateau(1.0f, 20) }; + + const float result{ math::PercentOvershoot(v, 0.0f) }; + + EXPECT_NEAR(result, 0.0f, math::Tolerance()); +} + +TEST_F(TestStepResponseMetrics, peak_time_finds_maximum_index) +{ + Vec v; + for (std::size_t i = 0; i < N; ++i) + v.at(i, 0) = 1.0f; + v.at(15, 0) = 1.5f; + + const float result{ math::PeakTime(v) }; + + EXPECT_NEAR(result, 15.0f, math::Tolerance()); +} + +TEST_F(TestStepResponseMetrics, peak_time_with_dt_scales_to_real_time) +{ + Vec v; + for (std::size_t i = 0; i < N; ++i) + v.at(i, 0) = 1.0f; + v.at(15, 0) = 1.5f; + const float dt{ 0.005f }; + + const float result{ math::PeakTime(v, dt) }; + + EXPECT_NEAR(result, 15.0f * dt, math::Tolerance()); +} + +TEST_F(TestStepResponseMetrics, steady_state_error_zero_for_perfect_step) +{ + Vec v; + for (std::size_t i = 0; i < N; ++i) + v.at(i, 0) = 1.0f; + + const float result{ math::SteadyStateError(v, 1.0f) }; + + EXPECT_NEAR(result, 0.0f, 1e-4f); +} + +TEST_F(TestStepResponseMetrics, steady_state_error_known_offset) +{ + Vec v; + for (std::size_t i = 0; i < N; ++i) + v.at(i, 0) = 0.9f; + + const float result{ math::SteadyStateError(v, 1.0f) }; + + EXPECT_NEAR(result, 0.1f, 1e-4f); +} From 18283fd46929bf0b8e5a7dc6e1e4844c4012450c Mon Sep 17 00:00:00 2001 From: gfs Date: Wed, 29 Jul 2026 17:31:56 +0200 Subject: [PATCH 2/3] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- doc/analysis/Decibels.md | 12 ++++++------ doc/analysis/README.md | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/analysis/Decibels.md b/doc/analysis/Decibels.md index 08b23722..adcd77b4 100644 --- a/doc/analysis/Decibels.md +++ b/doc/analysis/Decibels.md @@ -43,12 +43,12 @@ Both are simple differences in decibel space, exploiting the logarithm identity ## Complexity Analysis -| Operation | Time | Space | Notes | -|-----------------|----------|-------|------------------------------------| -| `ToDecibels` | O(1) | O(1) | One `log10` + one `max` + one `mul` | -| `FromDecibels` | O(1) | O(1) | One `pow` | -| `AttenuationDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction | -| `RippleDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction | +| Operation | Time | Space | Notes | +|-----------------|------|-------|-------------------------------------| +| `ToDecibels` | O(1) | O(1) | One `log10` + one `max` + one `mul` | +| `FromDecibels` | O(1) | O(1) | One `pow` | +| `AttenuationDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction | +| `RippleDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction | No state, no buffers. All operations are pure functions. diff --git a/doc/analysis/README.md b/doc/analysis/README.md index f55564aa..9e7606c7 100644 --- a/doc/analysis/README.md +++ b/doc/analysis/README.md @@ -12,7 +12,7 @@ Signal analysis algorithms for frequency-domain decomposition and spectral estim | [Discrete Cosine Transform](DiscreteCosineTransform.md) | Real-valued frequency decomposition via cosine basis functions, computed through FFT | | [Signal Detectors](SignalDetectors.md) | Peak hold, zero-crossing counter, and RMS envelope detectors for real-time signal monitoring | | [Decibels](Decibels.md) | `ToDecibels` / `FromDecibels` conversion helpers with zero-floor guard, plus attenuation and ripple utilities | -| [Goertzel Algorithm](GoertzelAlgorithm.md) | Single-bin DFT via a second-order recurrence for O(N) tone detection with O(1) memory | +| [Goertzel Algorithm](GoertzelAlgorithm.md) | Single-bin DFT via a second-order recurrence for O(N) tone detection with O(1) memory | ## Sub-domains From b7f334575c44f48b95693e7cdd90cf601e7cd726 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Wed, 29 Jul 2026 15:43:53 +0000 Subject: [PATCH 3/3] fix sonar findings --- numerical/analysis/CMakeLists.txt | 1 - numerical/analysis/Decibels.cpp | 9 --------- numerical/analysis/Decibels.hpp | 7 ------- numerical/math/CMakeLists.txt | 1 - numerical/math/StepResponseMetrics.cpp | 10 ---------- numerical/math/StepResponseMetrics.hpp | 14 +++----------- 6 files changed, 3 insertions(+), 39 deletions(-) delete mode 100644 numerical/analysis/Decibels.cpp delete mode 100644 numerical/math/StepResponseMetrics.cpp diff --git a/numerical/analysis/CMakeLists.txt b/numerical/analysis/CMakeLists.txt index 13348ade..540ff3fb 100644 --- a/numerical/analysis/CMakeLists.txt +++ b/numerical/analysis/CMakeLists.txt @@ -25,7 +25,6 @@ target_sources(numerical.analysis PRIVATE numerical_add_coverage_sources(numerical.analysis ConvolutionCorrelation.cpp - Decibels.cpp FastFourierTransformRadix2Impl.cpp GoertzelAlgorithm.cpp PowerDensitySpectrum.cpp diff --git a/numerical/analysis/Decibels.cpp b/numerical/analysis/Decibels.cpp deleted file mode 100644 index 78cff0bd..00000000 --- a/numerical/analysis/Decibels.cpp +++ /dev/null @@ -1,9 +0,0 @@ -#include "numerical/analysis/Decibels.hpp" - -namespace analysis -{ - template float ToDecibels(float); - template float FromDecibels(float); - template float AttenuationDb(float, float); - template float RippleDb(float, float); -} diff --git a/numerical/analysis/Decibels.hpp b/numerical/analysis/Decibels.hpp index 86cdea98..f53b4056 100644 --- a/numerical/analysis/Decibels.hpp +++ b/numerical/analysis/Decibels.hpp @@ -47,11 +47,4 @@ namespace analysis static_assert(std::is_floating_point_v, "RippleDb supports floating-point types only"); return ToDecibels(maxRatio) - ToDecibels(minRatio); } - -#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD - extern template float ToDecibels(float); - extern template float FromDecibels(float); - extern template float AttenuationDb(float, float); - extern template float RippleDb(float, float); -#endif } diff --git a/numerical/math/CMakeLists.txt b/numerical/math/CMakeLists.txt index 9b10d554..4eaca3d1 100644 --- a/numerical/math/CMakeLists.txt +++ b/numerical/math/CMakeLists.txt @@ -35,7 +35,6 @@ numerical_add_coverage_sources(numerical.math Matrix.cpp QNumber.cpp Quaternion.cpp - StepResponseMetrics.cpp ) add_subdirectory(test) diff --git a/numerical/math/StepResponseMetrics.cpp b/numerical/math/StepResponseMetrics.cpp deleted file mode 100644 index 8ac9f69a..00000000 --- a/numerical/math/StepResponseMetrics.cpp +++ /dev/null @@ -1,10 +0,0 @@ -#include "numerical/math/StepResponseMetrics.hpp" - -namespace math -{ - template float RiseTime(const Vector&, float, float); - template float SettlingTime(const Vector&, float, float, float); - template float PercentOvershoot(const Vector&, float); - template float PeakTime(const Vector&, float); - template float SteadyStateError(const Vector&, float); -} diff --git a/numerical/math/StepResponseMetrics.hpp b/numerical/math/StepResponseMetrics.hpp index 3cc0b0ec..1091ed09 100644 --- a/numerical/math/StepResponseMetrics.hpp +++ b/numerical/math/StepResponseMetrics.hpp @@ -18,8 +18,8 @@ namespace math { static_assert(std::is_floating_point_v, "RiseTime supports floating-point types only"); - const T lo{ T{ 0.1 } * steady }; - const T hi{ T{ 0.9 } * steady }; + const T lo{ T{ 0.1f } * steady }; + const T hi{ T{ 0.9f } * steady }; std::size_t iLo{ 0 }; bool foundLo{ false }; @@ -40,7 +40,7 @@ namespace math } template - [[nodiscard]] OPTIMIZE_FOR_SPEED T SettlingTime(const Vector& response, T steady, T band = T{ 0.02 }, T dt = T{ 1 }) + [[nodiscard]] OPTIMIZE_FOR_SPEED T SettlingTime(const Vector& response, T steady, T band = T{ 0.02f }, T dt = T{ 1 }) { static_assert(std::is_floating_point_v, "SettlingTime supports floating-point types only"); @@ -113,12 +113,4 @@ namespace math const T tailMean{ sum / static_cast(TailSize) }; return reference - tailMean; } - -#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD - extern template float RiseTime(const Vector&, float, float); - extern template float SettlingTime(const Vector&, float, float, float); - extern template float PercentOvershoot(const Vector&, float); - extern template float PeakTime(const Vector&, float); - extern template float SteadyStateError(const Vector&, float); -#endif }