feat: add lms adaptive filter

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                                        |
 | [Controllers](doc/controllers/README.md)                           | Bang-Bang/Hysteresis, PID, LQR, MPC, Saturation, Rate Limiter, Slew-Limited Saturation, Feedforward/2-DOF, Gain-Scheduled Controller, Lead-Lag Compensator |
 | [Dynamics](doc/dynamics/README.md)                                 | Euler-Lagrange, Newton-Euler, Recursive Newton-Euler, ABA            |
-| [Estimators](doc/estimators/README.md)                             | Linear Regression, Polynomial Fitting, Yule-Walker (offline), Recursive Least Squares (online) |
+| [Estimators](doc/estimators/README.md)                             | Linear Regression, Polynomial Fitting, Yule-Walker (offline), Recursive Least Squares, LMS / NLMS Adaptive Filter (online) |
 | [Filters](doc/filters/README.md)                                   | Kalman, Extended Kalman, Unscented Kalman, Alpha-Beta/Alpha-Beta-Gamma, FIR, IIR, Exponential Moving Average, Moving Average, Complementary, Median Filter, CIC (Cascaded Integrator-Comb), Notch/Comb Filter |
 | [Kinematics](doc/kinematics/README.md)                             | Forward Kinematics                                                   |
 | [Neural Network](doc/neural_network/README.md)                     | Layers, activations, losses, model                                   |

ROADMAP.md

@@ -30,7 +30,6 @@ Difficulty legend:
 | 15 | Biquad / Second-Order-Section cascade                | `filters/passive`         | ★★★☆☆      |
 | 19 | Luenberger observer + pole placement (Ackermann)     | `controllers`             | ★★★☆☆      |
 | 20 | Integral / servo state feedback (LQI)                | `controllers`             | ★★★☆☆      |
-| 21 | LMS / NLMS adaptive filter                           | `estimators/online`       | ★★★☆☆      |
 | 22 | Savitzky-Golay filter                                | `filters/passive`         | ★★★☆☆      |
 | 23 | CORDIC                                               | `math`                    | ★★★☆☆      |
 | 24 | Runge-Kutta ODE integrators (RK4 + Dormand-Prince)   | `solvers`                 | ★★★☆☆      |
@@ -166,12 +165,6 @@ Difficulty legend:
 - **Algorithm / paper:** B. D. O. Anderson, J. B. Moore, *Optimal Control: Linear Quadratic Methods* (1990).
 - **Reuses:** [Lqr.hpp](numerical/controllers/implementations/Lqr.hpp), [DARE](numerical/solvers/DiscreteAlgebraicRiccatiEquation.hpp).
 
-### 21. LMS / NLMS adaptive filter
-- **What:** Least-Mean-Squares (and normalized) adaptive FIR with steepest-descent weight update.
-- **Embedded value:** Adaptive noise cancellation, echo cancellation, system identification, active vibration control — a DSP staple.
-- **Algorithm / paper:** B. Widrow, M. Hoff, "Adaptive switching circuits," *IRE WESCON*, 1960; S. Haykin, *Adaptive Filter Theory*.
-- **Reuses:** [RecursiveLeastSquares.hpp](numerical/estimators/online/RecursiveLeastSquares.hpp) patterns, `math::Vector`.
-
 ### 22. Savitzky-Golay filter
 - **What:** Convolution smoother that fits a local polynomial, preserving peak height/width and giving smoothed derivatives.
 - **Embedded value:** Spectroscopy, ECG/PPG, and any signal where peaks matter and moving-average distortion is unacceptable.

doc/estimators/LmsAdaptiveFilter.md

@@ -0,0 +1,100 @@
+# LMS / NLMS Adaptive Filter
+
+## Overview & Motivation
+
+Real-time systems frequently need filters whose characteristics are unknown in advance or change over time. A fixed-coefficient FIR filter cannot cancel acoustic echo from an unknown loudspeaker-to-microphone path, equalize a time-varying channel, or identify the impulse response of a vibrating structure it has never seen before. The Least Mean Squares (LMS) algorithm and its normalized variant (NLMS) address this by adjusting a FIR filter's coefficients online, one sample at a time, without requiring matrix inversions or covariance storage.
+
+## Mathematical Theory
+
+### Prerequisites
+
+Let $\mathbf{w} \in \mathbb{R}^{N}$ be the weight vector (FIR coefficients) and $\mathbf{x}_n \in \mathbb{R}^{N}$ the tapped-delay-line vector at time $n$, with $x_n$ as the newest sample and $x_{n-N+1}$ the oldest. The filter output is
+
+$$y_n = \mathbf{w}_n^\top \mathbf{x}_n.$$
+
+Given a desired signal $d_n$, the estimation error is
+
+$$e_n = d_n - y_n.$$
+
+### Core Definitions
+
+The instantaneous squared error $J_n = e_n^2$ is a quadratic function of $\mathbf{w}_n$. Its gradient with respect to $\mathbf{w}$ is
+
+$$\nabla J_n = -2 e_n \mathbf{x}_n.$$
+
+### Derivation
+
+LMS takes a noisy steepest-descent step using the instantaneous gradient:
+
+$$\mathbf{w}_{n+1} = \mathbf{w}_n - \frac{\mu}{2} \nabla J_n = \mathbf{w}_n + \mu e_n \mathbf{x}_n.$$
+
+NLMS divides the step by the input energy to make the effective step size input-level independent:
+
+$$\mathbf{w}_{n+1} = \mathbf{w}_n + \frac{\mu}{\epsilon + \|\mathbf{x}_n\|^2} e_n \mathbf{x}_n,$$
+
+where $\epsilon > 0$ is a small regularizer that prevents division by zero during silence ($\|\mathbf{x}_n\|^2 \approx 0$).
+
+### Proof of Correctness
+
+Under the independence assumption (successive input vectors are uncorrelated), the expected weight vector converges to the Wiener solution $\mathbf{w}^* = \mathbf{R}^{-1}\mathbf{p}$ (where $\mathbf{R}$ is the input autocorrelation matrix and $\mathbf{p}$ is the cross-correlation vector) provided
+
+$$0 < \mu < \frac{2}{\lambda_{\max}(\mathbf{R})},$$
+
+a sufficient condition being $\mu < 2 / (N \cdot \sigma_x^2)$. For NLMS the stability range simplifies to $0 < \mu < 2$.
+
+## Complexity Analysis
+
+| Case | Time   | Space  | Notes                                       |
+|------|--------|--------|---------------------------------------------|
+| Any  | $O(N)$ | $O(N)$ | Two dot products + one AXPY, all length $N$ |
+
+No covariance matrix is stored; memory is exactly two vectors of length $N$ (weights and delay line).
+
+## Step-by-Step Walkthrough
+
+Consider a 2-tap system with true plant $\mathbf{w}^* = [0.5,\ 0.5]^\top$, $\mu = 0.1$, and a white input sequence.
+
+| Step | Input | Delay line $\mathbf{x}$ | Output $y$ | Desired $d$ | Error $e$ | Weight update direction     |
+|------|-------|-------------------------|------------|-------------|-----------|-----------------------------|
+| 1    | 1.0   | [1.0, 0.0]              | 0.0        | 0.5         | 0.5       | $+0.05\mathbf{x}$           |
+| 2    | 0.8   | [0.8, 1.0]              | 0.04       | 0.9         | 0.86      | $+0.086\mathbf{x}$          |
+| …    | …     | …                       | …          | …           | …         | converges to $\mathbf{w}^*$ |
+
+After many iterations the error approaches zero and the weights stabilize near $[0.5,\ 0.5]^\top$.
+
+## Pitfalls & Edge Cases
+
+- **Stability bound:** exceeding $\mu_{\max}$ causes exponential weight growth. For stationary white input, $\mu_{\max} = 2 / (N\sigma_x^2)$.
+- **Misadjustment:** even at steady state, random gradient noise keeps the weights oscillating around $\mathbf{w}^*$, adding excess mean-squared error proportional to $\mu N$.
+- **Slow convergence vs stability tradeoff:** small $\mu$ gives low misadjustment but slow tracking; large $\mu$ tracks fast but overshoots.
+- **Silent input:** $\|\mathbf{x}_n\|^2 = 0$ causes a singularity in plain NLMS; $\epsilon$ prevents this but freezes adaptation during silence.
+- **Non-persistent excitation:** if the input does not excite all directions, some weights may drift. Leaky-LMS ($\mathbf{w} \leftarrow (1-\rho)\mathbf{w} + \mu e \mathbf{x}$) bounds the drift.
+
+## Variants & Generalizations
+
+- **NLMS:** normalizes step by input energy; removes input-power dependence; recommended for signals with variable amplitude.
+- **Leaky LMS:** adds a forgetting term $\rho \ll 1$ to bound weight norm under insufficient excitation.
+- **Block LMS:** accumulates gradients over a block before updating; lower update rate but more FFT-friendly for long filters.
+- **Sign-LMS / Sign-Error LMS:** replaces $e$ or $\mathbf{x}$ with their signs; extremely simple hardware implementation at the cost of convergence rate.
+- **RLS:** replaces the stochastic gradient with the exact least-squares update; converges in $O(N)$ steps but requires $O(N^2)$ memory and time per update.
+
+## Applications
+
+- **Acoustic echo cancellation:** adapt the loudspeaker-to-microphone path in teleconferencing.
+- **Adaptive noise cancellation:** estimate and subtract a correlated noise source (e.g. engine noise on an aircraft headset).
+- **Channel equalization:** track a time-varying ISI channel in wireless modems.
+- **Active noise control / vibration control:** generate an anti-phase signal to cancel mechanical vibration.
+- **Online system identification:** recover the impulse response of an unknown plant in real time.
+
+## Connections to Other Algorithms
+
+- **RecursiveLeastSquares:** the exact-gradient counterpart; converges in far fewer samples but needs $O(N^2)$ storage and time per step. Prefer RLS for short filter lengths or fast convergence; prefer LMS/NLMS for embedded, resource-constrained paths.
+- **FIR (fixed-coefficient):** LMS turns a fixed FIR into a self-tuning one; the underlying convolution is identical.
+- **ConvolutionCorrelation:** used to compute the cross-correlation $\mathbf{p}$ in the batch Wiener solution; LMS approximates this online.
+
+## References & Further Reading
+
+- B. Widrow, M. E. Hoff, "Adaptive switching circuits," *IRE WESCON Convention Record*, 1960.
+- S. Haykin, *Adaptive Filter Theory*, 5th ed., Pearson, 2014 — chapters 5–8.
+- S. Sayed, *Fundamentals of Adaptive Filtering*, Wiley, 2003.
+- D. G. Manolakis, V. K. Ingle, S. M. Kogon, *Statistical and Adaptive Signal Processing*, Artech House, 2005.

numerical/estimators/online/CMakeLists.txt

@@ -12,10 +12,12 @@ target_link_libraries(numerical.estimators.online ${NUMERICAL_VISIBILITY}
 
 target_sources(numerical.estimators.online PRIVATE
     RecursiveLeastSquares.hpp
+    LmsAdaptiveFilter.hpp
 )
 
 numerical_add_coverage_sources(numerical.estimators.online
     RecursiveLeastSquares.cpp
+    LmsAdaptiveFilter.cpp
 )
 
 add_subdirectory(test)

numerical/estimators/online/LmsAdaptiveFilter.cpp

@@ -0,0 +1,6 @@
+#include "numerical/estimators/online/LmsAdaptiveFilter.hpp"
+
+namespace estimators
+{
+    template class LmsAdaptiveFilter<float, 4>;
+}

numerical/estimators/online/LmsAdaptiveFilter.hpp

@@ -0,0 +1,122 @@
+#pragma once
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC optimize("O3", "fast-math")
+#endif
+
+#include "numerical/math/CompilerOptimizations.hpp"
+#include <array>
+#include <cstddef>
+#include <type_traits>
+
+namespace estimators
+{
+    template<typename T, std::size_t Taps>
+    class LmsAdaptiveFilter
+    {
+        static_assert(std::is_floating_point_v<T>, "LmsAdaptiveFilter supports floating-point types only");
+        static_assert(Taps > 0, "Taps must be greater than zero");
+
+    public:
+        struct Result
+        {
+            T output{};
+            T error{};
+        };
+
+        explicit LmsAdaptiveFilter(T mu, bool normalized = false, T epsilon = T{ 1e-6 });
+
+        OPTIMIZE_FOR_SPEED Result Update(T input, T desired);
+        const std::array<T, Taps>& Weights() const;
+        void Reset();
+
+    private:
+        std::array<T, Taps> weights{};
+        std::array<T, Taps> delayLine{};
+        std::size_t head{ 0 };
+        T stepSize;
+        bool useNlms;
+        T regularizer;
+
+        T DotProduct() const;
+        T InputEnergy() const;
+        T GetSample(std::size_t tapIndex) const;
+    };
+
+    template<typename T, std::size_t Taps>
+    LmsAdaptiveFilter<T, Taps>::LmsAdaptiveFilter(T mu, bool normalized, T epsilon)
+        : stepSize{ mu }
+        , useNlms{ normalized }
+        , regularizer{ epsilon }
+    {
+    }
+
+    template<typename T, std::size_t Taps>
+    T LmsAdaptiveFilter<T, Taps>::GetSample(std::size_t tapIndex) const
+    {
+        std::size_t index{ (head + Taps - 1u - tapIndex) % Taps };
+        return delayLine[index];
+    }
+
+    template<typename T, std::size_t Taps>
+    T LmsAdaptiveFilter<T, Taps>::DotProduct() const
+    {
+        T result{ T{ 0 } };
+        for (std::size_t i{ 0 }; i < Taps; ++i)
+            result += weights[i] * GetSample(i);
+        return result;
+    }
+
+    template<typename T, std::size_t Taps>
+    T LmsAdaptiveFilter<T, Taps>::InputEnergy() const
+    {
+        T result{ T{ 0 } };
+        for (std::size_t i{ 0 }; i < Taps; ++i)
+        {
+            T s{ delayLine[i] };
+            result += s * s;
+        }
+        return result;
+    }
+
+    template<typename T, std::size_t Taps>
+    OPTIMIZE_FOR_SPEED typename LmsAdaptiveFilter<T, Taps>::Result LmsAdaptiveFilter<T, Taps>::Update(T input, T desired)
+    {
+        delayLine[head] = input;
+        head = (head + 1) % Taps;
+
+        T output{ DotProduct() };
+        T error{ desired - output };
+
+        T step{ stepSize };
+        if (useNlms)
+            step = stepSize / (regularizer + InputEnergy());
+
+        T scaled{ step * error };
+        for (std::size_t i{ 0 }; i < Taps; ++i)
+            weights[i] += scaled * GetSample(i);
+
+        return Result{ output, error };
+    }
+
+    template<typename T, std::size_t Taps>
+    const std::array<T, Taps>& LmsAdaptiveFilter<T, Taps>::Weights() const
+    {
+        return weights;
+    }
+
+    template<typename T, std::size_t Taps>
+    void LmsAdaptiveFilter<T, Taps>::Reset()
+    {
+        weights.fill(T{ 0 });
+        delayLine.fill(T{ 0 });
+        head = 0;
+    }
+}
+
+#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
+namespace estimators
+{
+    extern template class LmsAdaptiveFilter<float, 4>;
+}
+#endif

numerical/estimators/online/test/CMakeLists.txt

@@ -9,4 +9,5 @@ target_link_libraries(numerical.estimators.online_test PUBLIC
 
 target_sources(numerical.estimators.online_test PRIVATE
     TestRecursiveLeastSquares.cpp
+    TestLmsAdaptiveFilter.cpp
 )