Skip to content

Commit 59dc321

Browse files
Merge branch 'main' into feature/add-matrix-norms
2 parents 15c3581 + b61f468 commit 59dc321

14 files changed

Lines changed: 607 additions & 25 deletions

README.md

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal
1616

1717
| Category | Description |
1818
|--------------------------------------------------------------------|----------------------------------------------------------------------|
19-
| [Analysis](doc/analysis/README.md) | FFT, Real-Input FFT (RFFT), Power Spectral Density, DCT, Window Functions, Signal Detectors, Convolution & Correlation, Goertzel Algorithm |
19+
| [Analysis](doc/analysis/README.md) | FFT, Real-Input FFT (RFFT), Power Spectral Density, DCT, Window Functions, Signal Detectors, Convolution & Correlation, Goertzel Algorithm, Decibels |
2020
| [Control Analysis](doc/control_analysis/README.md) | Frequency Response, Root Locus, Controllability/Observability Matrices & Gramians |
2121
| [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 |
2222
| [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
2626
| [Neural Network](doc/neural_network/README.md) | Layers, activations, losses, model |
2727
| [Optimization](doc/optimization/README.md) | Gradient Descent |
2828
| [Regularization](doc/regularization/README.md) | L1 (Lasso), L2 (Ridge) |
29-
| [Math](doc/math/README.md) | CORDIC, Quaternion, MatrixNorms |
29+
| [Math](doc/math/README.md) | CORDIC, Quaternion, MatrixNorms, Step Response Metrics |
3030
| [Solvers](doc/solvers/README.md) | Gaussian Elimination, Levinson-Durbin, Durand-Kerner, Cholesky, DARE, Runge-Kutta ODE Integrators (RK4 + Dormand-Prince), Condition Number |
3131
| [Performance Optimization](doc/performance-optimization/README.md) | Compiler optimizations, SIMD |
3232

ROADMAP.md

Lines changed: 0 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,6 @@ Difficulty legend:
4848
| 45 | IIR filter design (Butterworth/Chebyshev + bilinear) | `filters/passive` | ★★★★★ |
4949
| 46 | H∞ state-feedback control | `robust_control` (new) | ★★★★★ |
5050
| 47 | Model Reference Adaptive Control (MRAC) | `nonlinear_control` (new) | ★★★★★ |
51-
| 48 | Decibel & magnitude-response helpers | `analysis` | ★☆☆☆☆ |
52-
| 49 | Step / transient-response metrics | `math` | ★★☆☆☆ |
5351
| 51 | Spectral radius / discrete stability margin | `math` | ★★★☆☆ |
5452
| 52 | Estimator consistency metrics (NEES / NIS) | `estimators` | ★★★☆☆ |
5553

@@ -560,23 +558,6 @@ and [`control_analysis/FrequencyResponse`](numerical/control_analysis/FrequencyR
560558
magnitude/phase. The items below are the missing pieces. All are **float-only**, no-heap, and operate
561559
on bounded `math::Vector`/`math::Matrix` inputs; tests are `TEST_F` on `float`.
562560

563-
### 48. Decibel & magnitude-response helpers ★☆☆☆☆ — `analysis`
564-
- **What:** `ToDecibels(ratio)` / `FromDecibels`, plus convenience magnitude(-in-dB) and attenuation
565-
helpers layered over `control_analysis::FrequencyResponse`.
566-
- **Metric value:** M2 (frequency response) — pass-band ripple / stop-band attenuation in dB, the
567-
natural unit for filter and controller tests.
568-
- **Algorithm:** `20·log10(·)`; guard the zero/`-inf` case with a floor.
569-
- **Reuses:** `math::TrigonometricFunctions`/`std::log10`, existing `FrequencyResponse`.
570-
571-
### 49. Step / transient-response metrics ★★☆☆☆ — `math`
572-
- **What:** From a bounded step-response `Vector`: `RiseTime` (10–90 %), `SettlingTime` (±band),
573-
`PercentOvershoot`, `PeakTime`, `SteadyStateError`.
574-
- **Metric value:** M3 (time-response) — the core acceptance criteria for every controller
575-
(`controllers/`), `LinearTimeInvariant`, and IIR filter.
576-
- **Algorithm:** single forward pass over the sampled response against the reference/steady value;
577-
standard control-systems definitions.
578-
- **Reuses:** `math::Vector`, `math::Statistics` for the steady-state estimate.
579-
580561
### 51. Spectral radius / discrete stability margin ★★★☆☆ — `math`
581562
- **What:** Dominant `|eigenvalue|` of a square (state/companion) matrix; `IsSchurStable` (all
582563
`|λ| < 1`) and the stability margin `1 − ρ(A)`.

doc/analysis/Decibels.md

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
# Decibels & Magnitude-Response Helpers
2+
3+
## Overview & Motivation
4+
5+
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.
6+
7+
## Mathematical Theory
8+
9+
### Magnitude Decibel Conversion
10+
11+
For a positive amplitude ratio $r > 0$, the equivalent level in decibels is:
12+
13+
$$L_{\mathrm{dB}} = 20 \log_{10}(r)$$
14+
15+
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.
16+
17+
### Inverse Conversion
18+
19+
$$r = 10^{L_{\mathrm{dB}}/20}$$
20+
21+
This inverse is exact for all finite $L_{\mathrm{dB}}$; no guard is needed on the output side.
22+
23+
### Zero and Negative Input Guard
24+
25+
$\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:
26+
27+
$$L_{\mathrm{dB}} = \max\!\left(20\log_{10}(r),\; L_{\min}\right), \quad r > 0$$
28+
$$L_{\mathrm{dB}} = L_{\min}, \quad r \leq 0$$
29+
30+
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.
31+
32+
### Derived Helpers
33+
34+
**Stop-band attenuation** between a pass-band ratio $r_p$ and a stop-band ratio $r_s$:
35+
36+
$$A = L_{\mathrm{dB}}(r_p) - L_{\mathrm{dB}}(r_s)$$
37+
38+
**Pass-band ripple** between the maximum and minimum in-band ratios $r_{\max}$ and $r_{\min}$:
39+
40+
$$\Delta = L_{\mathrm{dB}}(r_{\max}) - L_{\mathrm{dB}}(r_{\min})$$
41+
42+
Both are simple differences in decibel space, exploiting the logarithm identity $\log(a/b) = \log a - \log b$.
43+
44+
## Complexity Analysis
45+
46+
| Operation | Time | Space | Notes |
47+
|-----------------|------|-------|-------------------------------------|
48+
| `ToDecibels` | O(1) | O(1) | One `log10` + one `max` + one `mul` |
49+
| `FromDecibels` | O(1) | O(1) | One `pow` |
50+
| `AttenuationDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction |
51+
| `RippleDb` | O(1) | O(1) | Two `ToDecibels` + one subtraction |
52+
53+
No state, no buffers. All operations are pure functions.
54+
55+
## Step-by-Step Walkthrough
56+
57+
Converting a ratio of 10 to decibels:
58+
59+
1. Input $r = 10$; guard passes ($r > 0$).
60+
2. Compute $20 \cdot \log_{10}(10) = 20 \cdot 1 = 20$.
61+
3. Apply floor: $\max(20, -160) = 20$.
62+
4. Output: $20\,\mathrm{dB}$.
63+
64+
Round-trip for $r = 0.5$:
65+
66+
1. `ToDecibels(0.5)` = $20 \cdot \log_{10}(0.5) \approx -6.0206\,\mathrm{dB}$.
67+
2. `FromDecibels(-6.0206)` = $10^{-6.0206/20} \approx 0.5$.
68+
69+
## Pitfalls & Edge Cases
70+
71+
- Passing $r = 0$ produces $-\infty$ from `log10`; the floor guard prevents propagation into downstream computations.
72+
- 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.
73+
- 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.
74+
- `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.
75+
76+
## Variants & Generalizations
77+
78+
- 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.
79+
- 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.
80+
- 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.
81+
82+
## Applications
83+
84+
- Filter specification: stop-band attenuation and pass-band ripple in dB are the primary acceptance criteria for IIR/FIR designs.
85+
- Frequency response plots: `FrequencyResponse::Calculate()` already returns magnitude in dB using $20 \log_{10}$; these helpers provide the same conversion for ad-hoc analysis.
86+
- Controller gain margin is expressed in dB; converting from a linear ratio with `ToDecibels` avoids duplication.
87+
- Audio dynamic processing (compressor thresholds, limiter ceilings) and acoustic measurement both use dB natively.
88+
89+
## Connections to Other Algorithms
90+
91+
- `control_analysis::FrequencyResponse` internally applies $20 \log_{10}(\|H\|)$ on its magnitude output vector; these helpers are the scalar equivalent exposed for library consumers.
92+
- Pass-band ripple computed by `RippleDb` feeds directly into filter-design acceptance testing alongside the step/transient-response metrics in the evaluation primitives family.
93+
94+
## References & Further Reading
95+
96+
- Proakis, J. & Manolakis, D., "Digital Signal Processing", 4th ed., Prentice Hall, 2007 — Appendix A (decibel notation).
97+
- Zolzer, U., "DAFX: Digital Audio Effects", 2nd ed., Wiley, 2011 — Chapter 2 (level and gain in dB).
98+
- IEC 61672-1:2013, "Electroacoustics — Sound level meters — Part 1: Specifications."

doc/analysis/README.md

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,8 @@ Signal analysis algorithms for frequency-domain decomposition and spectral estim
1111
| [Power Spectral Density](PowerDensitySpectrum.md) | Estimation of signal power distribution across frequencies using Welch's method |
1212
| [Discrete Cosine Transform](DiscreteCosineTransform.md) | Real-valued frequency decomposition via cosine basis functions, computed through FFT |
1313
| [Signal Detectors](SignalDetectors.md) | Peak hold, zero-crossing counter, and RMS envelope detectors for real-time signal monitoring |
14-
| [Goertzel Algorithm](GoertzelAlgorithm.md) | Single-bin DFT via a second-order recurrence for O(N) tone detection with O(1) memory |
14+
| [Decibels](Decibels.md) | `ToDecibels` / `FromDecibels` conversion helpers with zero-floor guard, plus attenuation and ripple utilities |
15+
| [Goertzel Algorithm](GoertzelAlgorithm.md) | Single-bin DFT via a second-order recurrence for O(N) tone detection with O(1) memory |
1516

1617
## Sub-domains
1718

doc/math/README.md

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ Core mathematical primitives for numerical computation.
66

77
| Algorithm | Description |
88
|-----------------------------|-----------------------------------------------------------------------------------------------|
9-
| [CORDIC](Cordic.md) | Iterative shift-add engine for sin/cos, atan2, magnitude, and vector rotation — no multiplier |
10-
| [Quaternion](Quaternion.md) | Unit-quaternion rotation type: Hamilton product, SLERP, rotation-matrix and Euler conversions |
11-
| [MatrixNorms](MatrixNorms.md) | Frobenius, 1-norm, infinity-norm on matrices; vector L2 norm/normalize |
9+
| [CORDIC](Cordic.md) | Iterative shift-add engine for sin/cos, atan2, magnitude, and vector rotation — no multiplier |
10+
| [Quaternion](Quaternion.md) | Unit-quaternion rotation type: Hamilton product, SLERP, rotation-matrix and Euler conversions |
11+
| [MatrixNorms](MatrixNorms.md) | Frobenius, 1-norm, infinity-norm on matrices; vector L2 norm/normalize |
12+
| [Step Response Metrics](StepResponseMetrics.md) | Rise time, settling time, percent overshoot, peak time, and steady-state error from a bounded step-response vector |

doc/math/StepResponseMetrics.md

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
# Step / Transient-Response Metrics
2+
3+
## Overview & Motivation
4+
5+
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.
6+
7+
## Mathematical Theory
8+
9+
### Definitions
10+
11+
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$.
12+
13+
**Rise Time** $T_r$
14+
15+
The elapsed time for the response to travel from 10 % to 90 % of steady state:
16+
17+
$$T_r = (k_{90} - k_{10})\,\Delta t$$
18+
19+
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}\}$.
20+
21+
**Settling Time** $T_s$
22+
23+
The first time after which the response remains permanently inside the band $[(1-\delta)y_{ss},\,(1+\delta)y_{ss}]$ (typically $\delta = 0.02$):
24+
25+
$$T_s = (k^* + 1)\,\Delta t, \quad k^* = \max\{k : |y[k] - y_{ss}| > \delta\,|y_{ss}|\}$$
26+
27+
**Percent Overshoot** $\%OS$
28+
29+
$$\%OS = 100\,\frac{y_{\max} - y_{ss}}{y_{ss}}, \quad y_{\max} = \max_k y[k]$$
30+
31+
For an underdamped second-order system with damping ratio $\zeta$:
32+
33+
$$\%OS = 100\,\exp\!\left(-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}\right)$$
34+
35+
**Peak Time** $T_p$
36+
37+
$$T_p = k_p\,\Delta t, \quad k_p = \arg\max_k y[k]$$
38+
39+
For a continuous underdamped second-order system with natural frequency $\omega_n$:
40+
41+
$$T_p = \frac{\pi}{\omega_n\sqrt{1-\zeta^2}}$$
42+
43+
**Steady-State Error** $e_{ss}$
44+
45+
$$e_{ss} = r - \bar{y}_{\text{tail}}$$
46+
47+
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.
48+
49+
## Complexity Analysis
50+
51+
| Case | Time | Space | Notes |
52+
|---------|----------|--------|--------------------------------------------|
53+
| All | $O(N)$ | $O(1)$ | Single forward pass; no auxiliary storage |
54+
55+
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.
56+
57+
## Step-by-Step Walkthrough
58+
59+
Consider a 10-sample ramp to $y_{ss} = 1$ followed by a constant plateau (N = 20):
60+
61+
```
62+
k: 0 1 2 3 4 5 6 7 8 9 10 11 …
63+
y: 0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1 1 …
64+
```
65+
66+
- **Rise Time:** $k_{10} = 1$ (first sample $\ge 0.1$), $k_{90} = 9$ (first sample $\ge 0.9$). $T_r = 8\,\Delta t$.
67+
- **Settling Time:** With $\delta = 0.02$, last sample outside the band is $k = 9$. $T_s = 10\,\Delta t$.
68+
- **Percent Overshoot:** $y_{\max} = 1.0 = y_{ss}$, so $\%OS = 0$.
69+
- **Peak Time:** $k_p = 10$ (first occurrence of max). $T_p = 10\,\Delta t$.
70+
- **Steady-State Error:** Tail mean $= 1.0$, reference $= 1.0$. $e_{ss} = 0$.
71+
72+
## Pitfalls & Edge Cases
73+
74+
**Zero steady state.** Division by $y_{ss}$ in percent overshoot is guarded; the function returns zero when $y_{ss} = 0$ to avoid a NaN.
75+
76+
**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.
77+
78+
**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.
79+
80+
**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.
81+
82+
**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.
83+
84+
## Variants & Generalizations
85+
86+
- **Delay Time** $T_d$: the time to reach 50 % of steady state — obtainable with the same threshold-scan pattern.
87+
- **Band-relative rise time**: using a band other than 10–90 % (e.g., 20–80 %) is a trivial parameter change.
88+
- **Multi-channel:** applying the scalar functions element-wise to each row of a response matrix generalises to MIMO systems without algorithmic change.
89+
90+
## Applications
91+
92+
- Automated controller tuning acceptance: verify that a PID or LQR design meets specification ($T_r < T_{r,\text{spec}}$, $\%OS < \%OS_{\text{spec}}$, etc.).
93+
- Filter characterisation: measure the transient of a step fed through an IIR or FIR filter.
94+
- Hardware-in-the-loop test harnesses: compute metrics directly from sampled actuator responses.
95+
96+
## Connections to Other Algorithms
97+
98+
- **Statistics** (this library): the tail-mean for steady-state error replicates the `Mean` function on a sub-range.
99+
- **LinearTimeInvariant**: the primary source of step responses whose metrics are evaluated here.
100+
- **Filters/active** (Kalman, EKF): step-excitation tests use these metrics to validate estimator transient behaviour.
101+
- **Controllers**: PID and LQR tuning loops iterate until all five metrics satisfy design targets.
102+
103+
## References & Further Reading
104+
105+
- K. J. Åström and R. M. Murray, *Feedback Systems: An Introduction for Scientists and Engineers*, Princeton University Press, 2008. Chapter 10.
106+
- G. F. Franklin, J. D. Powell, and A. Emami-Naeini, *Feedback Control of Dynamic Systems*, 8th ed., Pearson, 2019. Chapter 3.
107+
- N. S. Nise, *Control Systems Engineering*, 8th ed., Wiley, 2019. Chapter 4.

numerical/analysis/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ target_link_libraries(numerical.analysis ${NUMERICAL_VISIBILITY}
1313

1414
target_sources(numerical.analysis PRIVATE
1515
ConvolutionCorrelation.hpp
16+
Decibels.hpp
1617
DiscreteCosineTransform.hpp
1718
FastFourierTransform.hpp
1819
FastFourierTransformRadix2Impl.hpp

0 commit comments

Comments
 (0)