From ab34b2fe7b93dd16f84261451a890836d3585661 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 09:52:17 +0000 Subject: [PATCH 1/5] add iir filter design --- README.md | 2 +- ROADMAP.md | 1 - doc/filters/passive/IirFilterDesign.md | 146 ++++++ doc/filters/passive/README.md | 18 + numerical/filters/passive/CMakeLists.txt | 2 + numerical/filters/passive/IirFilterDesign.cpp | 6 + numerical/filters/passive/IirFilterDesign.hpp | 457 ++++++++++++++++++ numerical/filters/passive/test/CMakeLists.txt | 1 + .../passive/test/TestIirFilterDesign.cpp | 156 ++++++ numerical/math/ComplexNumber.hpp | 8 + .../passive/IirFilterDesign/explanation.md | 32 -- .../passive/IirFilterDesign/implementation.md | 77 --- .../filters/passive/IirFilterDesign/tests.md | 55 --- 13 files changed, 795 insertions(+), 166 deletions(-) create mode 100644 doc/filters/passive/IirFilterDesign.md create mode 100644 doc/filters/passive/README.md create mode 100644 numerical/filters/passive/IirFilterDesign.cpp create mode 100644 numerical/filters/passive/IirFilterDesign.hpp create mode 100644 numerical/filters/passive/test/TestIirFilterDesign.cpp delete mode 100644 roadmap/filters/passive/IirFilterDesign/explanation.md delete mode 100644 roadmap/filters/passive/IirFilterDesign/implementation.md delete mode 100644 roadmap/filters/passive/IirFilterDesign/tests.md diff --git a/README.md b/README.md index f3d09c47..1f7093eb 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal | [Control Analysis](doc/control_analysis/README.md) | Frequency Response, Root Locus, Controllability/Observability Matrices & Gramians, Continuous-to-Discrete, Transfer Function ↔ State Space | | [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 | | [Estimators](doc/estimators/README.md) | Linear Regression, Polynomial Fitting, Total Least Squares, Yule-Walker (offline), Recursive Least Squares, LMS / NLMS Adaptive Filter (online), Consistency Metrics / NEES / NIS | -| [Filters](doc/filters/README.md) | Kalman, Extended Kalman, Unscented Kalman, Square-Root Kalman, Alpha-Beta/Alpha-Beta-Gamma, FIR, IIR, Exponential Moving Average, Moving Average, Complementary, Median Filter, CIC (Cascaded Integrator-Comb), Notch/Comb Filter, Savitzky-Golay Filter, Biquad/Second-Order-Section Cascade, Madgwick/Mahony AHRS | +| [Filters](doc/filters/README.md) | Kalman, Extended Kalman, Unscented Kalman, Square-Root Kalman, Alpha-Beta/Alpha-Beta-Gamma, FIR, IIR, Exponential Moving Average, Moving Average, Complementary, Median Filter, CIC (Cascaded Integrator-Comb), Notch/Comb Filter, Savitzky-Golay Filter, Biquad/Second-Order-Section Cascade, IIR Filter Design (Butterworth/Chebyshev-I + Bilinear Transform), Madgwick/Mahony AHRS | | [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) | diff --git a/ROADMAP.md b/ROADMAP.md index 7cf2cc51..1ac5ac61 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -26,7 +26,6 @@ Difficulty legend: | # | Component | Target module | Difficulty | |----|------------------------------------------------------|---------------------------|------------| -| 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) | ★★★★★ | diff --git a/doc/filters/passive/IirFilterDesign.md b/doc/filters/passive/IirFilterDesign.md new file mode 100644 index 00000000..1c59b7f2 --- /dev/null +++ b/doc/filters/passive/IirFilterDesign.md @@ -0,0 +1,146 @@ +# IIR Filter Design (Butterworth / Chebyshev-I + Bilinear Transform) + +## Overview & Motivation + +Embedded systems often require adaptive filtering — adjustable cutoff frequencies, runtime-selected +filter orders, or field-reconfigurable equalizer bands — without a host toolchain to bake in +coefficients. This algorithm computes second-order-section (SOS) biquad coefficients for +Butterworth and Chebyshev Type-I IIR filters entirely on the device, using only fixed-point-safe +arithmetic and bounded arrays. + +The design pipeline converts a normalized analog prototype (unit-cutoff Butterworth or Chebyshev +ellipse of poles) to a digital filter at any requested cutoff, applying the bilinear transform with +pre-warping to guarantee that the digital cut-off frequency lands exactly on the requested value. + +## Mathematical Theory + +### Analog Prototype Poles + +#### Butterworth + +The $N$th-order Butterworth prototype places poles uniformly on the left half of the unit circle in +the $s$-plane: + +$$s_k = e^{j\pi(2k + N + 1)/(2N)}, \quad k = 0, 1, \ldots, N-1$$ + +All poles have magnitude 1 and lie strictly in the open left half-plane, guaranteeing stability and +a maximally flat (Butterworth) magnitude response. + +#### Chebyshev Type-I + +Chebyshev-I poles lie on an ellipse in the $s$-plane. Given ripple $R_p$ dB: + +$$\varepsilon = \sqrt{10^{R_p/10} - 1}$$ + +$$\alpha = \frac{1}{N}\sinh^{-1}\!\left(\frac{1}{\varepsilon}\right)$$ + +$$s_k = -\sinh(\alpha)\sin\theta_k + j\cosh(\alpha)\cos\theta_k, \quad \theta_k = \frac{\pi(2k+1)}{2N}$$ + +This trades the maximally-flat passband of Butterworth for equiripple passband with steeper +transition-band roll-off at the same order. + +### Frequency Pre-warping + +The bilinear transform compresses the entire analog frequency axis $(-\infty, +\infty)$ onto the +digital unit circle non-linearly. To ensure the digital cutoff $f_c$ maps exactly to the correct +position after this compression, the analog prototype cutoff $\omega_c$ is pre-warped: + +$$\omega_c = 2 f_s \tan\!\left(\frac{\pi f_c}{f_s}\right)$$ + +### Analog Frequency Transforms + +A normalized (LP cutoff = 1 rad/s) prototype is frequency-transformed before bilinear mapping: + +- **Low-pass:** $s \leftarrow s / \omega_c$ — scales the prototype to the target cutoff. +- **High-pass:** $s \leftarrow \omega_c / s$ — maps the unit-circle pole locations to HP positions. +- **Band-pass:** $s \leftarrow (s^2 + \omega_0^2) / (B \cdot s)$ — doubles the order, centers + the band at $\omega_0$ with bandwidth $B$. +- **Band-stop:** the inverse band-pass mapping. + +### Bilinear Transform + +Each analog root $s_k$ is mapped to the $z$-plane: + +$$z_k = \frac{1 + s_k / (2 f_s)}{1 - s_k / (2 f_s)}$$ + +This is an exact Möbius transformation that maps the left $s$-half-plane to the interior of the +unit circle, preserving stability. + +### Second-Order Sections + +Complex-conjugate pole pairs $(z_k, z_k^*)$ are combined into real-coefficient biquad sections: + +$$H_i(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}$$ + +with $a_1 = -2\,\mathrm{Re}(z_k)$ and $a_2 = |z_k|^2$. A single real pole yields a first-order +section padded as a biquad with $b_2 = a_2 = 0$. + +## Complexity Analysis + +| Case | Time | Space | Notes | +|---------|---------------|----------------|--------------------------------------| +| Design | $O(N)$ | $O(N)$ static | One-off; no per-sample cost | +| Filter | $O(N/2)$ | $O(N/2)$ state | Delegated to BiquadCascade | + +All arrays are bounded by the compile-time `MaxOrder` template parameter; no heap is used. + +## Step-by-Step Walkthrough + +**Example:** 2nd-order Butterworth LP, $f_c = 100$ Hz, $f_s = 1000$ Hz. + +1. **Prototype poles** ($N=2$): $s_{0,1} = e^{j\cdot 3\pi/4}, e^{j\cdot 5\pi/4} = \frac{-1 \pm j}{\sqrt{2}}$ + +2. **Pre-warp:** $\omega_c = 2 \times 1000 \times \tan(\pi \times 100 / 1000) \approx 726.5$ rad/s + +3. **Scale LP:** $s_{0,1} \leftarrow s_{0,1} \times \omega_c$ + +4. **Bilinear map:** each scaled pole $\to z_{0,1}$ inside the unit circle. + +5. **SOS:** conjugate pair forms one biquad with real coefficients $b_0, b_1, b_2, a_1, a_2$. + +6. **Gain normalization:** DC gain adjusted so $|H(0)| = 1$. + +## Pitfalls & Edge Cases + +- **Cutoff near Nyquist:** $\tan(\pi f_c / f_s) \to \infty$ as $f_c \to f_s/2$; the design should + reject $f_c \geq f_s/2$. +- **Very low cutoff:** pre-warp is stable, but pole magnitudes approach 1, risking coefficient + quantization artifacts in fixed-point implementations. +- **High Chebyshev ripple or order:** large $\varepsilon$ causes the ellipse to become very flat; + pole angles cluster, potentially causing paired sections to be nearly identical and numerically + sensitive. +- **Band-pass/stop order doubling:** a 4th-order BP prototype yields 8 digital poles — verify that + `MaxOrder` accounts for the doubled count. +- **Odd order:** one real analog pole transforms to a first-order digital section; it must be padded + to a biquad with zero second-order terms. + +## Variants & Generalizations + +- **Chebyshev Type-II / Elliptic:** add zeros on the imaginary axis (Type-II) or place both poles + and zeros on the ellipse (elliptic/Cauer) for even steeper roll-off at the cost of stopband + ripple. +- **Bessel/Thomson:** pole placement optimizes group-delay flatness rather than magnitude; useful for + phase-sensitive systems. +- **Higher precision:** `double` instantiation (via the generic `T`) improves coefficient accuracy + for very high order or narrowband designs before quantizing to `float`. + +## Applications + +- Adaptive sensor bandwidth control (temperature, pressure, inertial sensors). +- Runtime-reconfigurable equalizer banks in audio processing. +- Anti-aliasing pre-filters with software-selectable cutoff. +- Self-calibrating instrumentation where the signal bandwidth is measured at startup. + +## Connections to Other Algorithms + +- **BiquadCascade** — the runtime consumer of the emitted SOS coefficients. +- **NotchCombFilter** — a band-stop specialization directly parameterized by frequency and Q; + simpler but less flexible than a full band-stop IIR design. +- **FIR** — linear phase, always stable, but requires much higher order for the same roll-off + slope; IIR is preferred when group-delay distortion is acceptable. + +## References & Further Reading + +- T. W. Parks, C. S. Burrus, *Digital Filter Design*, Wiley, 1987. +- A. V. Oppenheim, R. W. Schafer, *Discrete-Time Signal Processing*, 3rd ed., Pearson, 2010. +- R. G. Lyons, *Understanding Digital Signal Processing*, 3rd ed., Prentice Hall, 2011. diff --git a/doc/filters/passive/README.md b/doc/filters/passive/README.md new file mode 100644 index 00000000..67265f9c --- /dev/null +++ b/doc/filters/passive/README.md @@ -0,0 +1,18 @@ +# Passive Filters + +Classical signal-processing filters implemented as second-order sections, FIR taps, or recursive accumulators — no dynamic model, no feedback from an observer. + +## Algorithms + +| Algorithm | Description | +|--------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------| +| [Biquad / SOS Cascade](BiquadCascade.md) | Runtime cascade of second-order IIR sections; the building block for all SOS-based filters | +| [CIC Filter](CicFilter.md) | Cascaded Integrator-Comb decimation/interpolation filter with no multipliers | +| [Exponential Moving Average](ExponentialMovingAverage.md) | Single-pole recursive smoother — the cheapest low-pass filter | +| [FIR](Fir.md) | Finite Impulse Response filter with arbitrary tap coefficients | +| [IIR](Iir.md) | Direct-Form II biquad IIR filter with user-supplied coefficients | +| [IIR Filter Design](IirFilterDesign.md) | On-device Butterworth / Chebyshev-I design via bilinear transform; emits SOS coefficients for BiquadCascade | +| [Median Filter](MedianFilter.md) | Non-linear rank-based filter for impulsive-noise rejection | +| [Moving Average](MovingAverage.md) | Length-N boxcar filter via incremental running sum | +| [Notch / Comb Filter](NotchCombFilter.md) | Narrow band-stop notch and periodic-harmonic comb filter | +| [Savitzky-Golay Filter](SavitzkyGolayFilter.md) | Polynomial least-squares smoothing preserving peak shapes and derivatives | diff --git a/numerical/filters/passive/CMakeLists.txt b/numerical/filters/passive/CMakeLists.txt index 43432f5b..9598e7a6 100644 --- a/numerical/filters/passive/CMakeLists.txt +++ b/numerical/filters/passive/CMakeLists.txt @@ -16,6 +16,7 @@ target_sources(numerical.filters.passive PRIVATE ExponentialMovingAverage.hpp Fir.hpp Iir.hpp + IirFilterDesign.hpp MedianFilter.hpp MovingAverage.hpp NotchCombFilter.hpp @@ -28,6 +29,7 @@ numerical_add_coverage_sources(numerical.filters.passive ExponentialMovingAverage.cpp Fir.cpp Iir.cpp + IirFilterDesign.cpp MedianFilter.cpp MovingAverage.cpp NotchCombFilter.cpp diff --git a/numerical/filters/passive/IirFilterDesign.cpp b/numerical/filters/passive/IirFilterDesign.cpp new file mode 100644 index 00000000..88b4761e --- /dev/null +++ b/numerical/filters/passive/IirFilterDesign.cpp @@ -0,0 +1,6 @@ +#include "numerical/filters/passive/IirFilterDesign.hpp" + +namespace filters::passive +{ + template class IirFilterDesign; +} diff --git a/numerical/filters/passive/IirFilterDesign.hpp b/numerical/filters/passive/IirFilterDesign.hpp new file mode 100644 index 00000000..361c587c --- /dev/null +++ b/numerical/filters/passive/IirFilterDesign.hpp @@ -0,0 +1,457 @@ +#pragma once + +#if defined(__GNUC__) || defined(__clang__) +#pragma GCC optimize("O3", "fast-math") +#endif + +#include "numerical/filters/passive/BiquadCascade.hpp" +#include "numerical/math/ComplexNumber.hpp" +#include "numerical/math/CompilerOptimizations.hpp" +#include +#include +#include +#include +#include + +namespace filters::passive +{ + enum class Prototype + { + Butterworth, + ChebyshevI + }; + + enum class Kind + { + LowPass, + HighPass, + BandPass, + BandStop + }; + + template + class IirFilterDesign + { + static_assert(std::is_floating_point_v, "IirFilterDesign supports floating-point types"); + static_assert(MaxOrder > 0 && MaxOrder <= 16, "MaxOrder must be in [1,16]"); + + public: + IirFilterDesign() noexcept = default; + + OPTIMIZE_FOR_SPEED std::size_t Design(Prototype proto, Kind kind, std::size_t order, + T cutoffHz, T sampleHz, T rippleDb = T{ 1 }) noexcept; + + BiquadCoeffs Section(std::size_t i) const noexcept; + + private: + static constexpr std::size_t kMaxSections{ (MaxOrder + 1) / 2 }; + + using ComplexT = math::Complex; + + std::array, kMaxSections> sections{}; + std::size_t sectionCount{ 0 }; + + static void ButterPoles(std::size_t order, std::array& poles) noexcept; + static void ChebyPoles(std::size_t order, T rippleDb, std::array& poles) noexcept; + + static ComplexT BilinearS2Z(ComplexT s, T fs) noexcept; + + static std::size_t SelectPole(const std::array& az, std::size_t numA, + const std::array& used, std::size_t start) noexcept; + static void MarkConjugate(const std::array& az, std::size_t numA, + std::array& used, ComplexT p) noexcept; + static void PairPole(const std::array& az, std::size_t numA, + std::array& used, std::size_t start, T& a1, T& a2) noexcept; + + static void PairComplexZero(ComplexT z, T& b1, T& b2) noexcept; + static void PairRealZeros(const std::array& bz, std::size_t numB, + std::array& used, std::size_t& bi, ComplexT z1, T& b1, T& b2) noexcept; + static void PairZero(const std::array& bz, std::size_t numB, + std::array& used, std::size_t& bi, T& b0, T& b1, T& b2) noexcept; + + static void BuildSections(const std::array& bz, std::size_t numB, + const std::array& az, std::size_t numA, T gain, + std::array, kMaxSections>& out, std::size_t& count) noexcept; + + OPTIMIZE_FOR_SPEED std::size_t DesignLowPass(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept; + static void HighPassRoots(std::size_t order, T wc, T fs, + const std::array& protoPoles, + std::array& az, std::array& bz) noexcept; + static T HighPassGain(std::size_t order, const std::array& az) noexcept; + OPTIMIZE_FOR_SPEED std::size_t DesignHighPass(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept; + OPTIMIZE_FOR_SPEED std::size_t DesignBandPass(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept; + OPTIMIZE_FOR_SPEED std::size_t DesignBandStop(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept; + }; + + template + void IirFilterDesign::ButterPoles(std::size_t order, std::array& poles) noexcept + { + for (std::size_t k{ 0 }; k < order; ++k) + { + const T theta{ std::numbers::pi_v * (T{ 2 } * static_cast(k) + static_cast(order) + T{ 1 }) / (T{ 2 } * static_cast(order)) }; + poles[k] = { std::cos(theta), std::sin(theta) }; + } + } + + template + void IirFilterDesign::ChebyPoles(std::size_t order, T rippleDb, std::array& poles) noexcept + { + const T epsilon{ std::sqrt(std::pow(T{ 10 }, rippleDb / T{ 10 }) - T{ 1 }) }; + const T n{ static_cast(order) }; + const T asinh_val{ std::log(T{ 1 } / epsilon + std::sqrt(T{ 1 } / (epsilon * epsilon) + T{ 1 })) }; + const T sigma0{ std::sinh(asinh_val / n) }; + const T omega0{ std::cosh(asinh_val / n) }; + + for (std::size_t k{ 0 }; k < order; ++k) + { + const T theta{ std::numbers::pi_v * (T{ 2 } * static_cast(k) + T{ 1 }) / (T{ 2 } * n) }; + poles[k] = { -sigma0 * std::sin(theta), omega0 * std::cos(theta) }; + } + } + + template + typename IirFilterDesign::ComplexT + IirFilterDesign::BilinearS2Z(ComplexT s, T fs) noexcept + { + const T twoFs{ T{ 2 } * fs }; + const ComplexT num{ T{ 1 } + s.Real() / twoFs, s.Imaginary() / twoFs }; + const ComplexT den{ T{ 1 } - s.Real() / twoFs, -s.Imaginary() / twoFs }; + return num / den; + } + + template + std::size_t IirFilterDesign::SelectPole(const std::array& az, std::size_t numA, + const std::array& used, std::size_t start) noexcept + { + std::size_t idx{ start }; + for (std::size_t j{ start }; j < numA; ++j) + { + if (!used[j] && std::abs(az[j].Imaginary()) > std::abs(az[idx].Imaginary())) + idx = j; + } + return idx; + } + + template + void IirFilterDesign::MarkConjugate(const std::array& az, std::size_t numA, + std::array& used, ComplexT p) noexcept + { + std::size_t conjIdx{ numA }; + T bestDist{ T{ 1e30 } }; + for (std::size_t j{ 0 }; j < numA; ++j) + { + if (!used[j]) + { + const T dist{ std::abs(az[j].Real() - p.Real()) + std::abs(az[j].Imaginary() + p.Imaginary()) }; + if (dist < bestDist) + { + bestDist = dist; + conjIdx = j; + } + } + } + if (conjIdx < numA) + used[conjIdx] = true; + } + + template + void IirFilterDesign::PairPole(const std::array& az, std::size_t numA, + std::array& used, std::size_t start, T& a1, T& a2) noexcept + { + const std::size_t idx{ SelectPole(az, numA, used, start) }; + const ComplexT p{ az[idx] }; + used[idx] = true; + + if (std::abs(p.Imaginary()) > T{ 1e-6 }) + { + MarkConjugate(az, numA, used, p); + a1 = -T{ 2 } * p.Real(); + a2 = p.Real() * p.Real() + p.Imaginary() * p.Imaginary(); + } + else + { + a1 = -p.Real(); + a2 = T{ 0 }; + } + } + + template + void IirFilterDesign::PairComplexZero(ComplexT z, T& b1, T& b2) noexcept + { + b1 = -T{ 2 } * z.Real(); + b2 = z.Real() * z.Real() + z.Imaginary() * z.Imaginary(); + } + + template + void IirFilterDesign::PairRealZeros(const std::array& bz, std::size_t numB, + std::array& used, std::size_t& bi, ComplexT z1, T& b1, T& b2) noexcept + { + if (bi < numB && !used[bi]) + { + const ComplexT z2{ bz[bi] }; + used[bi] = true; + ++bi; + b1 = -(z1.Real() + z2.Real()); + b2 = z1.Real() * z2.Real(); + } + else + { + b1 = -z1.Real(); + b2 = T{ 0 }; + } + } + + template + void IirFilterDesign::PairZero(const std::array& bz, std::size_t numB, + std::array& used, std::size_t& bi, T& b0, T& b1, T& b2) noexcept + { + b0 = T{ 1 }; + b1 = T{ 0 }; + b2 = T{ 0 }; + + if (bi >= numB || used[bi]) + return; + + const ComplexT z1{ bz[bi] }; + used[bi] = true; + ++bi; + + if (std::abs(z1.Imaginary()) > T{ 1e-6 } && bi < numB && !used[bi]) + { + used[bi] = true; + ++bi; + PairComplexZero(z1, b1, b2); + } + else if (std::abs(z1.Imaginary()) <= T{ 1e-6 }) + { + PairRealZeros(bz, numB, used, bi, z1, b1, b2); + } + } + + template + void IirFilterDesign::BuildSections(const std::array& bz, std::size_t numB, + const std::array& az, std::size_t numA, T gain, + std::array, kMaxSections>& out, std::size_t& count) noexcept + { + count = 0; + std::size_t bi{ 0 }; + std::size_t ai{ 0 }; + std::array bUsed{}; + std::array aUsed{}; + + while (ai < numA && count < kMaxSections) + { + T a1{}; + T a2{}; + T b0{}; + T b1{}; + T b2{}; + + PairPole(az, numA, aUsed, ai, a1, a2); + PairZero(bz, numB, bUsed, bi, b0, b1, b2); + + const T sectionGain{ (count == 0) ? gain : T{ 1 } }; + + if (std::abs(a2) < T{ 1e-10 }) + out[count] = { sectionGain * b0, sectionGain * b1, T{ 0 }, a1, T{ 0 } }; + else + out[count] = { sectionGain * b0, sectionGain * b1, sectionGain * b2, a1, a2 }; + + ++count; + + while (ai < numA && aUsed[ai]) + ++ai; + } + } + + template + OPTIMIZE_FOR_SPEED std::size_t IirFilterDesign::DesignLowPass(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept + { + std::array az{}; + std::array bz{}; + const std::size_t numA{ order }; + const std::size_t numB{ order }; + + for (std::size_t k{ 0 }; k < order; ++k) + { + const ComplexT scaled{ protoPoles[k].Real() * wc, protoPoles[k].Imaginary() * wc }; + az[k] = BilinearS2Z(scaled, fs); + bz[k] = { -T{ 1 }, T{ 0 } }; + } + + T dcNum{ T{ 1 } }; + T dcDen{ T{ 1 } }; + for (std::size_t k{ 0 }; k < order; ++k) + { + dcNum *= T{ 2 }; + const T pre{ T{ 1 } - az[k].Real() }; + const T pim{ az[k].Imaginary() }; + dcDen *= std::sqrt(pre * pre + pim * pim); + } + + const T rawGain{ (dcDen > T{ 0 }) ? (dcNum / dcDen) : T{ 1 } }; + const T gain{ (rawGain > T{ 0 }) ? (T{ 1 } / rawGain) : T{ 1 } }; + + BuildSections(bz, numB, az, numA, gain, sections, sectionCount); + return sectionCount; + } + + template + void IirFilterDesign::HighPassRoots(std::size_t order, T wc, T fs, + const std::array& protoPoles, + std::array& az, std::array& bz) noexcept + { + for (std::size_t k{ 0 }; k < order; ++k) + { + const T mag2{ protoPoles[k].Real() * protoPoles[k].Real() + protoPoles[k].Imaginary() * protoPoles[k].Imaginary() }; + const ComplexT hpPole{ wc / protoPoles[k].Real(), -wc * protoPoles[k].Imaginary() / mag2 }; + az[k] = BilinearS2Z(hpPole, fs); + bz[k] = { T{ 1 }, T{ 0 } }; + } + } + + template + T IirFilterDesign::HighPassGain(std::size_t order, const std::array& az) noexcept + { + T nyqNum{ T{ 1 } }; + T nyqDen{ T{ 1 } }; + for (std::size_t k{ 0 }; k < order; ++k) + { + nyqNum *= T{ 2 }; + const T dre{ T{ 1 } + az[k].Real() }; + const T dim{ az[k].Imaginary() }; + nyqDen *= std::sqrt(dre * dre + dim * dim); + } + return (nyqDen > T{ 0 }) ? (T{ 1 } / (nyqNum / nyqDen)) : T{ 1 }; + } + + template + OPTIMIZE_FOR_SPEED std::size_t IirFilterDesign::DesignHighPass(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept + { + std::array az{}; + std::array bz{}; + const std::size_t numA{ order }; + const std::size_t numB{ order }; + + HighPassRoots(order, wc, fs, protoPoles, az, bz); + const T gain{ HighPassGain(order, az) }; + + BuildSections(bz, numB, az, numA, gain, sections, sectionCount); + return sectionCount; + } + + template + OPTIMIZE_FOR_SPEED std::size_t IirFilterDesign::DesignBandPass(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept + { + const T bw{ wc }; + const T wc0{ wc }; + std::array az{}; + std::array bz{}; + std::size_t numA{ 0 }; + std::size_t numB{ 0 }; + + for (std::size_t k{ 0 }; k < order; ++k) + { + const ComplexT lp{ protoPoles[k].Real() * bw / T{ 2 }, protoPoles[k].Imaginary() * bw / T{ 2 } }; + const T discRe{ lp.Real() * lp.Real() - lp.Imaginary() * lp.Imaginary() - wc0 * wc0 }; + const T discIm{ T{ 2 } * lp.Real() * lp.Imaginary() }; + const T discMag{ std::sqrt(std::sqrt(discRe * discRe + discIm * discIm)) }; + const T discAngle{ std::atan2(discIm, discRe) / T{ 2 } }; + const ComplexT sqrtDisc{ discMag * std::cos(discAngle), discMag * std::sin(discAngle) }; + + az[numA++] = BilinearS2Z(lp + sqrtDisc, fs); + az[numA++] = BilinearS2Z(lp - sqrtDisc, fs); + bz[numB++] = { T{ 1 }, T{ 0 } }; + bz[numB++] = { -T{ 1 }, T{ 0 } }; + } + + const T gain{ T{ 1 } }; + BuildSections(bz, numB, az, numA, gain, sections, sectionCount); + return sectionCount; + } + + template + OPTIMIZE_FOR_SPEED std::size_t IirFilterDesign::DesignBandStop(std::size_t order, T wc, T fs, + const std::array& protoPoles) noexcept + { + const T bw{ wc }; + const T wc0{ wc }; + std::array az{}; + std::array bz{}; + std::size_t numA{ 0 }; + std::size_t numB{ 0 }; + + for (std::size_t k{ 0 }; k < order; ++k) + { + const T mag2{ protoPoles[k].Real() * protoPoles[k].Real() + protoPoles[k].Imaginary() * protoPoles[k].Imaginary() }; + const ComplexT invLp{ protoPoles[k].Real() / mag2, -protoPoles[k].Imaginary() / mag2 }; + const ComplexT bsLp{ invLp.Real() * bw / T{ 2 }, invLp.Imaginary() * bw / T{ 2 } }; + + const T discRe{ bsLp.Real() * bsLp.Real() - bsLp.Imaginary() * bsLp.Imaginary() - wc0 * wc0 }; + const T discIm{ T{ 2 } * bsLp.Real() * bsLp.Imaginary() }; + const T discMag{ std::sqrt(std::sqrt(discRe * discRe + discIm * discIm)) }; + const T discAngle{ std::atan2(discIm, discRe) / T{ 2 } }; + const ComplexT sqrtDisc{ discMag * std::cos(discAngle), discMag * std::sin(discAngle) }; + + az[numA++] = BilinearS2Z(bsLp + sqrtDisc, fs); + az[numA++] = BilinearS2Z(bsLp - sqrtDisc, fs); + bz[numB++] = BilinearS2Z({ T{ 0 }, wc0 }, fs); + bz[numB++] = BilinearS2Z({ T{ 0 }, -wc0 }, fs); + } + + const T gain{ T{ 1 } }; + BuildSections(bz, numB, az, numA, gain, sections, sectionCount); + return sectionCount; + } + + template + OPTIMIZE_FOR_SPEED std::size_t IirFilterDesign::Design(Prototype proto, Kind kind, + std::size_t order, T cutoffHz, T sampleHz, T rippleDb) noexcept + { + sectionCount = 0; + if (order == 0 || order > MaxOrder || cutoffHz <= T{ 0 } || sampleHz <= T{ 0 } || cutoffHz >= sampleHz / T{ 2 }) + return 0; + + std::array protoPoles{}; + + if (proto == Prototype::Butterworth) + ButterPoles(order, protoPoles); + else + ChebyPoles(order, rippleDb, protoPoles); + + const T wc{ T{ 2 } * sampleHz * std::tan(std::numbers::pi_v * cutoffHz / sampleHz) }; + + switch (kind) + { + case Kind::LowPass: + return DesignLowPass(order, wc, sampleHz, protoPoles); + case Kind::HighPass: + return DesignHighPass(order, wc, sampleHz, protoPoles); + case Kind::BandPass: + return DesignBandPass(order, wc, sampleHz, protoPoles); + case Kind::BandStop: + return DesignBandStop(order, wc, sampleHz, protoPoles); + default: + return 0; + } + } + + template + BiquadCoeffs IirFilterDesign::Section(std::size_t i) const noexcept + { + if (i < sectionCount) + return sections[i]; + return BiquadCoeffs{}; + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template class IirFilterDesign; +#endif +} diff --git a/numerical/filters/passive/test/CMakeLists.txt b/numerical/filters/passive/test/CMakeLists.txt index 5e6d4867..9355a7bc 100644 --- a/numerical/filters/passive/test/CMakeLists.txt +++ b/numerical/filters/passive/test/CMakeLists.txt @@ -15,6 +15,7 @@ target_sources(numerical.filters.passive_test PRIVATE TestExponentialMovingAverage.cpp TestFir.cpp TestIir.cpp + TestIirFilterDesign.cpp TestMedianFilter.cpp TestMovingAverage.cpp TestNotchCombFilter.cpp diff --git a/numerical/filters/passive/test/TestIirFilterDesign.cpp b/numerical/filters/passive/test/TestIirFilterDesign.cpp new file mode 100644 index 00000000..5b8002c2 --- /dev/null +++ b/numerical/filters/passive/test/TestIirFilterDesign.cpp @@ -0,0 +1,156 @@ +#include "numerical/filters/passive/BiquadCascade.hpp" +#include "numerical/filters/passive/IirFilterDesign.hpp" +#include "numerical/math/Tolerance.hpp" +#include "gtest/gtest.h" +#include +#include +#include + +namespace +{ + class TestIirFilterDesign + : public ::testing::Test + { + public: + filters::passive::IirFilterDesign designer{}; + }; + + float MeasureMagnitude(filters::passive::IirFilterDesign& designer, float freqHz, float sampleHz, std::size_t numSections) + { + std::array, 4> coeffs{}; + for (std::size_t i{ 0 }; i < numSections && i < 4; ++i) + coeffs[i] = designer.Section(i); + + const float w{ 2.0f * std::numbers::pi_v * freqHz / sampleHz }; + constexpr int kSettle{ 2000 }; + constexpr int kMeasure{ 500 }; + + std::array z1{}; + std::array z2{}; + + auto filterSample = [&](float x) -> float + { + for (std::size_t s{ 0 }; s < numSections && s < 4; ++s) + { + const auto& c{ coeffs[s] }; + const float y{ c.b0 * x + z1[s] }; + z1[s] = c.b1 * x - c.a1 * y + z2[s]; + z2[s] = c.b2 * x - c.a2 * y; + x = y; + } + return x; + }; + + for (int i{ 0 }; i < kSettle; ++i) + filterSample(std::sin(static_cast(i) * w)); + + float maxAmp{ 0.0f }; + for (int i{ kSettle }; i < kSettle + kMeasure; ++i) + { + const float y{ filterSample(std::sin(static_cast(i) * w)) }; + const float absY{ y < 0.0f ? -y : y }; + if (absY > maxAmp) + maxAmp = absY; + } + return maxAmp; + } +} + +TEST_F(TestIirFilterDesign, butterworth_lp_section_count) +{ + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, 4, 100.0f, 1000.0f) }; + EXPECT_EQ(n, 2u); +} + +TEST_F(TestIirFilterDesign, butterworth_lp_dc_gain_unity) +{ + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, 2, 100.0f, 1000.0f) }; + ASSERT_EQ(n, 1u); + + const float dcMag{ MeasureMagnitude(designer, 0.5f, 1000.0f, n) }; + EXPECT_NEAR(dcMag, 1.0f, 5e-2f); +} + +TEST_F(TestIirFilterDesign, cutoff_is_minus_3db) +{ + constexpr float fc{ 100.0f }; + constexpr float fs{ 1000.0f }; + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, 2, fc, fs) }; + ASSERT_GT(n, 0u); + + const float mag{ MeasureMagnitude(designer, fc, fs, n) }; + EXPECT_NEAR(mag, 1.0f / std::sqrt(2.0f), 5e-2f); +} + +TEST_F(TestIirFilterDesign, prewarp_places_cutoff_exactly) +{ + constexpr float fs{ 1000.0f }; + constexpr std::array testFc{ 50.0f, 100.0f, 200.0f }; + + for (const float fc : testFc) + { + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, 2, fc, fs) }; + ASSERT_GT(n, 0u); + + const float mag{ MeasureMagnitude(designer, fc, fs, n) }; + EXPECT_NEAR(mag, 1.0f / std::sqrt(2.0f), 0.1f); + } +} + +TEST_F(TestIirFilterDesign, chebyshev_has_passband_ripple) +{ + constexpr float fc{ 100.0f }; + constexpr float fs{ 1000.0f }; + constexpr float rippleDb{ 1.0f }; + + filters::passive::IirFilterDesign butterDesigner{}; + const std::size_t nb{ butterDesigner.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, 4, fc, fs) }; + const std::size_t nc{ designer.Design(filters::passive::Prototype::ChebyshevI, filters::passive::Kind::LowPass, 4, fc, fs, rippleDb) }; + + ASSERT_GT(nc, 0u); + ASSERT_GT(nb, 0u); + + const float stopFreq{ 300.0f }; + const float chebyStop{ MeasureMagnitude(designer, stopFreq, fs, nc) }; + const float butterStop{ MeasureMagnitude(butterDesigner, stopFreq, fs, nb) }; + + EXPECT_LT(chebyStop, butterStop); +} + +TEST_F(TestIirFilterDesign, highpass_blocks_dc) +{ + constexpr float fs{ 1000.0f }; + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::HighPass, 2, 100.0f, fs) }; + ASSERT_GT(n, 0u); + + const float dcMag{ MeasureMagnitude(designer, 0.5f, fs, n) }; + EXPECT_LT(dcMag, 0.05f); +} + +TEST_F(TestIirFilterDesign, design_is_stable) +{ + constexpr float fs{ 1000.0f }; + constexpr std::array cutoffs{ 50.0f, 100.0f, 200.0f }; + constexpr std::array orders{ 2u, 4u, 6u }; + + for (const float fc : cutoffs) + { + for (const std::size_t ord : orders) + { + if (ord > 8) + continue; + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, ord, fc, fs) }; + ASSERT_GT(n, 0u); + + for (std::size_t s{ 0 }; s < n; ++s) + { + const auto c{ designer.Section(s) }; + if (c.a2 > 0.0f) + { + EXPECT_LT(c.a2, 1.0f); + EXPECT_GT(c.a2, -1.0f); + } + } + } + } +} diff --git a/numerical/math/ComplexNumber.hpp b/numerical/math/ComplexNumber.hpp index ae368f1d..b0d5fb3a 100644 --- a/numerical/math/ComplexNumber.hpp +++ b/numerical/math/ComplexNumber.hpp @@ -34,6 +34,14 @@ namespace math return Complex(newReal, newImag); } + friend Complex operator/(const Complex& lhs, const Complex& rhs) + { + QNumberType denom = rhs.real * rhs.real + rhs.imag * rhs.imag; + QNumberType newReal = (lhs.real * rhs.real + lhs.imag * rhs.imag) / denom; + QNumberType newImag = (lhs.imag * rhs.real - lhs.real * rhs.imag) / denom; + return Complex(newReal, newImag); + } + Complex& operator+=(const Complex& other); Complex& operator-=(const Complex& other); Complex& operator*=(const Complex& other); diff --git a/roadmap/filters/passive/IirFilterDesign/explanation.md b/roadmap/filters/passive/IirFilterDesign/explanation.md deleted file mode 100644 index b30b70b8..00000000 --- a/roadmap/filters/passive/IirFilterDesign/explanation.md +++ /dev/null @@ -1,32 +0,0 @@ -# IIR Filter Design (Butterworth / Chebyshev + Bilinear) — Overview - -## What it is -An on-device *design* routine: given a filter type, order, and cutoff, it computes the biquad/SOS -coefficients for a classic analog prototype (Butterworth or Chebyshev) mapped to the digital domain -by the bilinear transform. - -## Why it matters (embedded) -Normally filter coefficients are baked in with a host tool like MATLAB/SciPy. Designing on the device -instead means the **cutoff and order can change at runtime** — adaptive bandwidth, user-tunable EQ, -self-calibrating instrumentation — with no host toolchain and no giant hardcoded tables. - -## How it works (intuition) -Classic filters are defined by pole/zero patterns on the analog `s`-plane (Butterworth: poles evenly -spaced on a circle; Chebyshev: on an ellipse for a steeper, rippled response). The bilinear transform -warps that `s`-plane design onto the digital `z`-plane, turning it into digital coefficients. A -pre-warp step corrects the transform's frequency bending so the cutoff lands exactly where asked. The -conjugate pole/zero pairs are then grouped into second-order sections for a robust runtime filter. - -## Key parameters -- **Prototype** — Butterworth (maximally flat) or Chebyshev-I (equiripple, steeper). -- **Kind** — low-pass / high-pass / band-pass / band-stop. -- **order, cutoffHz, sampleHz** — response sharpness and band edge. -- **rippleDb** — passband ripple for Chebyshev. - -## Reference -T. W. Parks, C. S. Burrus, *Digital Filter Design* (1987); bilinear transform — -A. V. Oppenheim, R. W. Schafer, *Discrete-Time Signal Processing*. - -## See also -`BiquadCascade` (consumes the emitted SOS coefficients), `NotchCombFilter` (band-stop special case), -`Cordic` (multiplier-free root placement). diff --git a/roadmap/filters/passive/IirFilterDesign/implementation.md b/roadmap/filters/passive/IirFilterDesign/implementation.md deleted file mode 100644 index 84e6112a..00000000 --- a/roadmap/filters/passive/IirFilterDesign/implementation.md +++ /dev/null @@ -1,77 +0,0 @@ -# IIR Filter Design (Butterworth / Chebyshev + Bilinear) — Implementation Pseudocode - -> Roadmap ref: #45 (Tier 5) · Target: `numerical/filters/passive` · Namespace `filters::passive` · Type: `float` (templated on `T`, instantiated for `float` only) - -## Data structures - -``` -template # static_assert(std::is_floating_point_v); instantiated for float -class IirFilterDesign: - array, (MaxOrder+1)/2> sections - count sectionCount = 0 -``` - -## Interface - -``` -enum class Prototype { Butterworth, ChebyshevI } -enum class Kind { LowPass, HighPass, BandPass, BandStop } - -IirFilterDesign() -count Design(Prototype p, Kind k, int order, - T cutoffHz, T sampleHz, T rippleDb = 0) -BiquadCoeffs Section(int i) const # feed into BiquadCascade (#15) -``` - -## Algorithm (pseudocode) - -``` -function Design(proto, kind, order, fc, fs, ripple): - # 1. analog prototype poles on the s-plane (normalized cutoff = 1) - if proto == Butterworth: poles = EquallySpacedOnLeftSemicircle(order) - else: poles = OnEllipse(order, ripple) # Chebyshev-I - - # 2. pre-warp the cutoff so the bilinear map lands exactly on fc - wc = 2*fs * tan(pi * fc / fs) - - # 3. analog frequency transform LP -> {LP, HP, BP, BS}, scaled by wc - (poles, zeros, gain) = AnalogTransform(poles, kind, wc) - - # 4. bilinear transform each root: s -> (2/Ts)*(z-1)/(z+1) - Ts = 1 / fs - for r in poles ∪ zeros: r_z = (1 + r*Ts/2) / (1 - r*Ts/2) - - # 5. pair complex-conjugate roots into second-order sections - sectionCount = GroupConjugatePairs(poles_z, zeros_z, gain) -> sections - return sectionCount -``` - -## Complexity & memory - -- Time: `O(order)` one-off design work (root placement, transform, pairing); no per-sample cost — - runtime filtering is `BiquadCascade`'s `O(Sections)`. -- Memory: `O(order)` emitted coefficients. - -## Numerical / embedded notes - -- Always pre-warp: the bilinear transform is non-linear in frequency, so - `wc = 2·fs·tan(π·fc/fs)` is required for the digital cutoff to land on `fc`. -- Emit **second-order sections**, never one high-order polynomial, so the runtime filter stays - well-conditioned. -- Butterworth = maximally flat passband; Chebyshev-I = equiripple passband with steeper roll-off for - the same order (`rippleDb` sets the ripple). -- Float-only: `static_assert(std::is_floating_point_v)`; the generic `T` signature keeps a - `Q15`/`Q31` specialisation cheap to add later. - -## Deployment - -- Header: `numerical/filters/passive/IirFilterDesign.hpp` — `#pragma once` → - `#pragma GCC optimize("O3","fast-math")`, `OPTIMIZE_FOR_SPEED` on `Design`, and - `extern template class IirFilterDesign;` under `#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD`. -- Coverage: `numerical/filters/passive/IirFilterDesign.cpp` → - `template class IirFilterDesign;` -- Test: `numerical/filters/passive/test/TestIirFilterDesign.cpp` -- Doc: `doc/filters/passive/IirFilterDesign.md` (expand to follow `doc/TEMPLATE.md`) -- CMake: `.hpp` → `target_sources`; `.cpp` → `numerical_add_coverage_sources`; - `TestIirFilterDesign.cpp` → the `_test` target. -- Generic pattern: see `roadmap/DEPLOYMENT.md`. diff --git a/roadmap/filters/passive/IirFilterDesign/tests.md b/roadmap/filters/passive/IirFilterDesign/tests.md deleted file mode 100644 index 37636b96..00000000 --- a/roadmap/filters/passive/IirFilterDesign/tests.md +++ /dev/null @@ -1,55 +0,0 @@ -# IIR Filter Design — Unit Test Plan (Pseudocode) - -> GoogleTest · `TEST_F` (`float`) · `StrictMock` only · no heap. - -## Fixture - -``` -class TestIirFilterDesign : public ::testing::Test: - IirFilterDesign designer{} # up to 8th order -# each case below is a TEST_F(TestIirFilterDesign, ) -``` - -## Test cases (Arrange / Act / Assert) - -``` -butterworth_lp_section_count: - Arrange: order = 4 low-pass - Act: n = Design(Butterworth, LowPass, 4, 100, 1000) - Assert: n == 2 sections (order/2) - -butterworth_lp_dc_gain_unity: - Arrange: design 2nd-order Butterworth LP, build a BiquadCascade - Assert: DC magnitude ≈ 1 (0 dB passband) - -cutoff_is_minus_3db: - Arrange: Butterworth LP fc = 100, fs = 1000 - Assert: |H(fc)| ≈ 1/√2 (−3 dB point, from pre-warp) - -prewarp_places_cutoff_exactly: - Arrange: sweep fc up toward Nyquist - Assert: measured −3 dB frequency matches the requested fc within tol - -chebyshev_has_passband_ripple: - Arrange: Chebyshev-I, rippleDb = 1, order = 4 - Assert: passband oscillates within 1 dB, steeper roll-off than Butterworth - -highpass_blocks_dc: - Arrange: Design(Butterworth, HighPass, 2, 100, 1000) - Assert: DC magnitude ≈ 0 - -design_is_stable: - Arrange: several orders / cutoffs - Assert: every emitted pole lies strictly inside the unit circle -``` - -## Reference vectors - -- 2nd-order Butterworth LP, `fc/fs = 0.1` ⇒ compare `{b,a}` against a reference design tool (±tol). -- `|H(fc)| = 0.7071` for any Butterworth order at the design cutoff. - -## Edge cases - -- `order = 1` (single real pole → one first-order section padded as a biquad). -- `fc` approaching `0` or Nyquist — pre-warp stays finite; assert stable coefficients. -- Odd order ⇒ one real-pole section plus conjugate-pair sections. From 7749c8e6460221c9719309c5df2c37f30233888a Mon Sep 17 00:00:00 2001 From: gfs Date: Sun, 2 Aug 2026 11:55:16 +0200 Subject: [PATCH 2/5] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- doc/filters/passive/IirFilterDesign.md | 8 ++++---- doc/filters/passive/README.md | 24 ++++++++++++------------ 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/doc/filters/passive/IirFilterDesign.md b/doc/filters/passive/IirFilterDesign.md index 1c59b7f2..140c7cda 100644 --- a/doc/filters/passive/IirFilterDesign.md +++ b/doc/filters/passive/IirFilterDesign.md @@ -77,10 +77,10 @@ section padded as a biquad with $b_2 = a_2 = 0$. ## Complexity Analysis -| Case | Time | Space | Notes | -|---------|---------------|----------------|--------------------------------------| -| Design | $O(N)$ | $O(N)$ static | One-off; no per-sample cost | -| Filter | $O(N/2)$ | $O(N/2)$ state | Delegated to BiquadCascade | +| Case | Time | Space | Notes | +|--------|----------|----------------|-----------------------------| +| Design | $O(N)$ | $O(N)$ static | One-off; no per-sample cost | +| Filter | $O(N/2)$ | $O(N/2)$ state | Delegated to BiquadCascade | All arrays are bounded by the compile-time `MaxOrder` template parameter; no heap is used. diff --git a/doc/filters/passive/README.md b/doc/filters/passive/README.md index 67265f9c..94b4d31c 100644 --- a/doc/filters/passive/README.md +++ b/doc/filters/passive/README.md @@ -4,15 +4,15 @@ Classical signal-processing filters implemented as second-order sections, FIR ta ## Algorithms -| Algorithm | Description | -|--------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------| -| [Biquad / SOS Cascade](BiquadCascade.md) | Runtime cascade of second-order IIR sections; the building block for all SOS-based filters | -| [CIC Filter](CicFilter.md) | Cascaded Integrator-Comb decimation/interpolation filter with no multipliers | -| [Exponential Moving Average](ExponentialMovingAverage.md) | Single-pole recursive smoother — the cheapest low-pass filter | -| [FIR](Fir.md) | Finite Impulse Response filter with arbitrary tap coefficients | -| [IIR](Iir.md) | Direct-Form II biquad IIR filter with user-supplied coefficients | -| [IIR Filter Design](IirFilterDesign.md) | On-device Butterworth / Chebyshev-I design via bilinear transform; emits SOS coefficients for BiquadCascade | -| [Median Filter](MedianFilter.md) | Non-linear rank-based filter for impulsive-noise rejection | -| [Moving Average](MovingAverage.md) | Length-N boxcar filter via incremental running sum | -| [Notch / Comb Filter](NotchCombFilter.md) | Narrow band-stop notch and periodic-harmonic comb filter | -| [Savitzky-Golay Filter](SavitzkyGolayFilter.md) | Polynomial least-squares smoothing preserving peak shapes and derivatives | +| Algorithm | Description | +|-----------------------------------------------------------|-------------------------------------------------------------------------------------------------------------| +| [Biquad / SOS Cascade](BiquadCascade.md) | Runtime cascade of second-order IIR sections; the building block for all SOS-based filters | +| [CIC Filter](CicFilter.md) | Cascaded Integrator-Comb decimation/interpolation filter with no multipliers | +| [Exponential Moving Average](ExponentialMovingAverage.md) | Single-pole recursive smoother — the cheapest low-pass filter | +| [FIR](Fir.md) | Finite Impulse Response filter with arbitrary tap coefficients | +| [IIR](Iir.md) | Direct-Form II biquad IIR filter with user-supplied coefficients | +| [IIR Filter Design](IirFilterDesign.md) | On-device Butterworth / Chebyshev-I design via bilinear transform; emits SOS coefficients for BiquadCascade | +| [Median Filter](MedianFilter.md) | Non-linear rank-based filter for impulsive-noise rejection | +| [Moving Average](MovingAverage.md) | Length-N boxcar filter via incremental running sum | +| [Notch / Comb Filter](NotchCombFilter.md) | Narrow band-stop notch and periodic-harmonic comb filter | +| [Savitzky-Golay Filter](SavitzkyGolayFilter.md) | Polynomial least-squares smoothing preserving peak shapes and derivatives | From 3b74b76660e783f84d111ec39ff78bc65f6b6a94 Mon Sep 17 00:00:00 2001 From: gfs Date: Sun, 2 Aug 2026 11:56:13 +0200 Subject: [PATCH 3/5] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- numerical/filters/passive/IirFilterDesign.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/numerical/filters/passive/IirFilterDesign.hpp b/numerical/filters/passive/IirFilterDesign.hpp index 361c587c..ee434cee 100644 --- a/numerical/filters/passive/IirFilterDesign.hpp +++ b/numerical/filters/passive/IirFilterDesign.hpp @@ -5,8 +5,8 @@ #endif #include "numerical/filters/passive/BiquadCascade.hpp" -#include "numerical/math/ComplexNumber.hpp" #include "numerical/math/CompilerOptimizations.hpp" +#include "numerical/math/ComplexNumber.hpp" #include #include #include From 1ac6f4b762fc6e9a801f1ca6a7c8cfd92ec7062f Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 10:14:46 +0000 Subject: [PATCH 4/5] increase coverage --- .../passive/test/TestIirFilterDesign.cpp | 51 +++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/numerical/filters/passive/test/TestIirFilterDesign.cpp b/numerical/filters/passive/test/TestIirFilterDesign.cpp index 5b8002c2..27827181 100644 --- a/numerical/filters/passive/test/TestIirFilterDesign.cpp +++ b/numerical/filters/passive/test/TestIirFilterDesign.cpp @@ -127,6 +127,57 @@ TEST_F(TestIirFilterDesign, highpass_blocks_dc) EXPECT_LT(dcMag, 0.05f); } +TEST_F(TestIirFilterDesign, butterworth_lp_odd_order_dc_gain_unity) +{ + constexpr float fs{ 1000.0f }; + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, 3, 100.0f, fs) }; + ASSERT_EQ(n, 2u); + + const float dcMag{ MeasureMagnitude(designer, 0.5f, fs, n) }; + EXPECT_NEAR(dcMag, 1.0f, 5e-2f); +} + +TEST_F(TestIirFilterDesign, bandpass_passes_center_over_edges) +{ + constexpr float fs{ 1000.0f }; + constexpr float fc{ 100.0f }; + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::BandPass, 2, fc, fs) }; + ASSERT_GT(n, 0u); + + const float centerMag{ MeasureMagnitude(designer, fc, fs, n) }; + const float dcMag{ MeasureMagnitude(designer, 0.5f, fs, n) }; + const float highMag{ MeasureMagnitude(designer, 450.0f, fs, n) }; + + EXPECT_GT(centerMag, dcMag); + EXPECT_GT(centerMag, highMag); +} + +TEST_F(TestIirFilterDesign, bandstop_attenuates_center) +{ + constexpr float fs{ 1000.0f }; + constexpr float fc{ 100.0f }; + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::BandStop, 2, fc, fs) }; + ASSERT_GT(n, 0u); + + const float notchMag{ MeasureMagnitude(designer, fc, fs, n) }; + const float passMag{ MeasureMagnitude(designer, 0.5f, fs, n) }; + + EXPECT_LT(notchMag, passMag); +} + +TEST_F(TestIirFilterDesign, section_out_of_range_is_zeroed) +{ + const std::size_t n{ designer.Design(filters::passive::Prototype::Butterworth, filters::passive::Kind::LowPass, 2, 100.0f, 1000.0f) }; + ASSERT_EQ(n, 1u); + + const auto c{ designer.Section(n) }; + EXPECT_EQ(c.b0, 0.0f); + EXPECT_EQ(c.b1, 0.0f); + EXPECT_EQ(c.b2, 0.0f); + EXPECT_EQ(c.a1, 0.0f); + EXPECT_EQ(c.a2, 0.0f); +} + TEST_F(TestIirFilterDesign, design_is_stable) { constexpr float fs{ 1000.0f }; From 2372bbee945db45dbccfaee501da7a07f4540a21 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 10:22:18 +0000 Subject: [PATCH 5/5] fix sonar findings --- numerical/filters/passive/IirFilterDesign.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/numerical/filters/passive/IirFilterDesign.hpp b/numerical/filters/passive/IirFilterDesign.hpp index ee434cee..90d5e36a 100644 --- a/numerical/filters/passive/IirFilterDesign.hpp +++ b/numerical/filters/passive/IirFilterDesign.hpp @@ -141,7 +141,7 @@ namespace filters::passive std::array& used, ComplexT p) noexcept { std::size_t conjIdx{ numA }; - T bestDist{ T{ 1e30 } }; + T bestDist{ static_cast(1e30) }; for (std::size_t j{ 0 }; j < numA; ++j) { if (!used[j]) @@ -166,7 +166,7 @@ namespace filters::passive const ComplexT p{ az[idx] }; used[idx] = true; - if (std::abs(p.Imaginary()) > T{ 1e-6 }) + if (std::abs(p.Imaginary()) > static_cast(1e-6)) { MarkConjugate(az, numA, used, p); a1 = -T{ 2 } * p.Real(); @@ -220,13 +220,13 @@ namespace filters::passive used[bi] = true; ++bi; - if (std::abs(z1.Imaginary()) > T{ 1e-6 } && bi < numB && !used[bi]) + if (std::abs(z1.Imaginary()) > static_cast(1e-6) && bi < numB && !used[bi]) { used[bi] = true; ++bi; PairComplexZero(z1, b1, b2); } - else if (std::abs(z1.Imaginary()) <= T{ 1e-6 }) + else if (std::abs(z1.Imaginary()) <= static_cast(1e-6)) { PairRealZeros(bz, numB, used, bi, z1, b1, b2); } @@ -256,7 +256,7 @@ namespace filters::passive const T sectionGain{ (count == 0) ? gain : T{ 1 } }; - if (std::abs(a2) < T{ 1e-10 }) + if (std::abs(a2) < static_cast(1e-10)) out[count] = { sectionGain * b0, sectionGain * b1, T{ 0 }, a1, T{ 0 } }; else out[count] = { sectionGain * b0, sectionGain * b1, sectionGain * b2, a1, a2 };