add hilbert transform

Author: gabrielfrasantos

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, Decibels |
+| [Analysis](doc/analysis/README.md)                                 | FFT, Real-Input FFT (RFFT), Power Spectral Density, DCT, Window Functions, Signal Detectors, Convolution & Correlation, Goertzel Algorithm, Decibels, Hilbert Transform / Analytic 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 |
 | [Dynamics](doc/dynamics/README.md)                                 | Euler-Lagrange, Newton-Euler, Recursive Newton-Euler, ABA            |

ROADMAP.md

@@ -29,7 +29,6 @@ Difficulty legend:
 |----|------------------------------------------------------|---------------------------|------------|
 | 35 | Disturbance Observer (DOB)                           | `robust_control` (new)    | ★★★★☆      |
 | 36 | Active Disturbance Rejection Control (ADRC + ESO)    | `robust_control` (new)    | ★★★★☆      |
-| 37 | Hilbert transform / analytic signal / envelope       | `analysis`                | ★★★★☆      |
 | 38 | Discrete Wavelet Transform (Haar / Daubechies)       | `analysis`                | ★★★★☆      |
 | 39 | Square-root / Information Kalman filter              | `filters/active`          | ★★★★☆      |
 | 40 | Feedback linearization                               | `nonlinear_control` (new) | ★★★★☆      |

doc/analysis/HilbertTransform.md

@@ -0,0 +1,103 @@
+# Hilbert Transform / Analytic Signal
+
+## Overview & Motivation
+
+Every real-valued signal contains both positive and negative frequency components that carry redundant information. The Hilbert transform discards the negative half and pairs the original signal with a version shifted by exactly −90° at every frequency, producing a **complex analytic signal** whose magnitude tracks instantaneous amplitude and whose angle tracks instantaneous phase.
+
+This is the standard technique for AM/DSB demodulation, vibration envelope detection, bearing-fault analysis, single-sideband modulation, and instantaneous-frequency measurement — all tasks where the underlying amplitude or phase varies slowly relative to a carrier and must be tracked in real time.
+
+## Mathematical Theory
+
+### Hilbert Transform
+
+For a real signal $x(t)$ the Hilbert transform is the convolution
+
+$$\mathcal{H}\{x\}(t) = \frac{1}{\pi} \, \text{p.v.} \int_{-\infty}^{\infty} \frac{x(\tau)}{t - \tau} \, d\tau$$
+
+which is equivalently a multiplication of each frequency component by $-j \operatorname{sgn}(f)$, i.e. a $-90°$ phase rotation for positive frequencies and $+90°$ for negative.
+
+### Analytic Signal
+
+The analytic signal is
+
+$$z(t) = x(t) + j\,\mathcal{H}\{x\}(t)$$
+
+Its one-sided spectrum satisfies $Z(f) = 0$ for $f < 0$, $Z(0) = X(0)$, and $Z(f) = 2X(f)$ for $f > 0$.
+
+### Instantaneous Attributes
+
+| Quantity | Formula |
+|----------|---------|
+| Amplitude (envelope) | $A(t) = |z(t)| = \sqrt{x^2 + \mathcal{H}^2\{x\}}$ |
+| Phase | $\phi(t) = \arg z(t) = \operatorname{atan2}(\mathcal{H}\{x\}, x)$ |
+| Frequency | $f_i(t) = \frac{1}{2\pi}\frac{d\phi}{dt}$ (phase unwrapped before differencing) |
+
+### FFT Method (Marple)
+
+Given the $N$-point DFT $X[k]$ of a real sequence, the analytic signal is recovered by:
+
+$$H[k] = \begin{cases} X[0] & k = 0 \\ 2X[k] & 1 \le k < N/2 \\ X[N/2] & k = N/2 \\ 0 & N/2 < k < N \end{cases}$$
+
+followed by the inverse DFT of $H$.
+
+### FIR Approximation
+
+A Type III/IV antisymmetric FIR filter with impulse response
+
+$$h[k] = \frac{2}{\pi k} w[k], \quad k \text{ odd}; \quad h[k] = 0, \quad k \text{ even}$$
+
+(with a Hamming window $w[k]$) approximates the ideal $-90°$ phase shift across its passband. The delayed original signal (center-tap copy) is paired with the filtered output to form the approximate analytic signal.
+
+## Complexity Analysis
+
+| Method | Time per call | Space | Notes |
+|--------|--------------|-------|-------|
+| FFT (block) | $O(N \log N)$ | $2N$ complex words | Exact, latency $N$ |
+| FIR (streaming) | $O(P)$ per sample | $P$ words state | Approx., latency $(P-1)/2$ |
+| Feature extraction | $O(1)$ | None | atan2, sqrt, subtract |
+
+$P$ = number of FIR taps.
+
+## Step-by-Step Walkthrough
+
+**FFT method on $N = 8$, input $x = \cos(2\pi n / 8)$:**
+
+1. Forward DFT yields $X[1] = 4$, $X[7] = 4$, all others zero.
+2. Apply one-sided weighting: $H[1] = 8$, $H[7] = 0$.
+3. Inverse DFT of $H$ produces the imaginary part $\sin(2\pi n / 8)$.
+4. Analytic signal: $z[n] = \cos(2\pi n/8) + j\sin(2\pi n/8)$, amplitude $\equiv 1$.
+
+## Pitfalls & Edge Cases
+
+- **DC and Nyquist bins** — their imaginary parts must remain zero; the one-sided formula preserves this.
+- **Phase unwrapping** — before computing instantaneous frequency, the phase difference must be mapped to $(-\pi, \pi]$ to suppress $2\pi$ jumps.
+- **Block-edge artefacts** — the FFT method treats the block as periodic; trim the first and last few samples when asserting accuracy.
+- **FIR group delay** — the FIR output is delayed by $(P-1)/2$ samples relative to the input; align before comparing to the FFT method.
+- **FIR bandwidth** — the FIR Hilbert approximation degrades near DC and Nyquist; use only over the filter's flat passband.
+
+## Variants & Generalizations
+
+- **Block FFT method** — exact, no ripple, requires a complete block; suitable for offline or buffered processing.
+- **FIR streaming method** — causal, constant memory, approximate; length and window choice trade accuracy against latency.
+- **Quadrature-oscillator method** — for narrowband signals, a simple IIR resonator can approximate the 90° shift with even lower cost.
+
+## Applications
+
+- AM/DSB envelope demodulation and AGC.
+- Vibration and bearing-fault analysis (envelope of resonance band).
+- Single-sideband (SSB) radio modulation/demodulation.
+- Instantaneous frequency tracking in FM receivers and Doppler radar.
+- Phase and frequency estimation in PLLs.
+
+## Connections to Other Algorithms
+
+- Uses `FastFourierTransform` / `FastFourierTransformRadix2Impl` as the block back-end.
+- `RealFastFourierTransform` is an alternative back-end for purely real inputs.
+- `ConvolutionCorrelation` provides the FIR convolution primitive used by the streaming path.
+- `SignalDetectors` (RMS envelope, peak hold) offers cheaper non-coherent envelope estimation.
+
+## References & Further Reading
+
+- S. L. Marple, "Computing the Discrete-Time Analytic Signal via FFT," *IEEE Transactions on Signal Processing*, 47(9), 2600–2603, 1999.
+- A. V. Oppenheim and R. W. Schafer, *Discrete-Time Signal Processing*, 3rd ed., Prentice Hall, Ch. 12.
+- S. W. Smith, *The Scientist and Engineer's Guide to Digital Signal Processing*, Ch. 9 (available free online).

doc/analysis/README.md

@@ -13,6 +13,7 @@ Signal analysis algorithms for frequency-domain decomposition and spectral estim
 | [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                         |
+| [Hilbert Transform](HilbertTransform.md)                | Analytic signal and instantaneous amplitude/phase/frequency via FFT one-sided spectrum or FIR approximation    |
 
 ## Sub-domains
 

numerical/analysis/CMakeLists.txt

@@ -18,6 +18,7 @@ target_sources(numerical.analysis PRIVATE
     FastFourierTransform.hpp
     FastFourierTransformRadix2Impl.hpp
     GoertzelAlgorithm.hpp
+    HilbertTransform.hpp
     PowerDensitySpectrum.hpp
     RealFastFourierTransform.hpp
     SignalDetectors.hpp
@@ -28,6 +29,7 @@ numerical_add_coverage_sources(numerical.analysis
     DiscreteCosineTransform.cpp
     FastFourierTransformRadix2Impl.cpp
     GoertzelAlgorithm.cpp
+    HilbertTransform.cpp
     PowerDensitySpectrum.cpp
     RealFastFourierTransform.cpp
     SignalDetectors.cpp

numerical/analysis/HilbertTransform.cpp

@@ -0,0 +1,7 @@
+#include "numerical/analysis/HilbertTransform.hpp"
+
+namespace analysis
+{
+    template class AnalyticSignalFft<float, 64>;
+    template class HilbertFir<float, 31>;
+}

numerical/analysis/HilbertTransform.hpp

@@ -0,0 +1,166 @@
+#pragma once
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC optimize("O3", "fast-math")
+#endif
+
+#include "infra/util/BoundedVector.hpp"
+#include "numerical/analysis/FastFourierTransform.hpp"
+#include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/ComplexNumber.hpp"
+#include <array>
+#include <cmath>
+#include <numbers>
+#include <type_traits>
+
+namespace analysis
+{
+    template<typename T, std::size_t N>
+    class AnalyticSignalFft
+    {
+        static_assert(std::is_floating_point_v<T>, "AnalyticSignalFft supports floating-point types");
+        static_assert(N >= 4 && (N & (N - 1)) == 0, "AnalyticSignalFft length N must be a power of two >= 4");
+
+    public:
+        using Complex = math::Complex<T>;
+        using VectorComplex = infra::BoundedVector<Complex>;
+        using VectorReal = infra::BoundedVector<T>;
+
+        explicit AnalyticSignalFft(FastFourierTransform<T>& fft);
+
+        OPTIMIZE_FOR_SPEED VectorComplex& Analytic(VectorReal& x);
+
+        static T InstantaneousAmplitude(Complex a);
+        static T InstantaneousPhase(Complex a);
+        static T InstantaneousFrequency(T phaseNow, T phasePrev, T ts);
+
+    private:
+        FastFourierTransform<T>& fft;
+        typename VectorComplex::template WithMaxSize<N> spectrum;
+        typename VectorComplex::template WithMaxSize<N> analytic;
+    };
+
+    template<typename T, std::size_t Taps>
+    class HilbertFir
+    {
+        static_assert(std::is_floating_point_v<T>, "HilbertFir supports floating-point types");
+        static_assert(Taps >= 3 && (Taps % 2) == 1, "HilbertFir requires an odd number of taps >= 3");
+
+    public:
+        using Complex = math::Complex<T>;
+
+        HilbertFir();
+
+        OPTIMIZE_FOR_SPEED Complex Filter(T x);
+
+    private:
+        static constexpr std::size_t centerTap{ (Taps - 1) / 2 };
+
+        std::array<T, Taps> coeff{};
+        std::array<T, Taps> delay{};
+        std::size_t writeIndex{ 0 };
+    };
+
+    /// AnalyticSignalFft implementation ///
+
+    template<typename T, std::size_t N>
+    AnalyticSignalFft<T, N>::AnalyticSignalFft(FastFourierTransform<T>& fft)
+        : fft{ fft }
+    {
+        spectrum.resize(N);
+        analytic.resize(N);
+    }
+
+    template<typename T, std::size_t N>
+    OPTIMIZE_FOR_SPEED typename AnalyticSignalFft<T, N>::VectorComplex& AnalyticSignalFft<T, N>::Analytic(VectorReal& x)
+    {
+        VectorComplex& X{ fft.Forward(x) };
+
+        spectrum[0] = X[0];
+        spectrum[N / 2] = X[N / 2];
+
+        for (std::size_t k{ 1 }; k < N / 2; ++k)
+            spectrum[k] = Complex{ T(2) * X[k].Real(), T(2) * X[k].Imaginary() };
+
+        for (std::size_t k{ N / 2 + 1 }; k < N; ++k)
+            spectrum[k] = Complex{ T(0), T(0) };
+
+        typename VectorComplex::template WithMaxSize<N> hilbertSpectrum{};
+        hilbertSpectrum.resize(N);
+        for (std::size_t k{ 0 }; k < N; ++k)
+            hilbertSpectrum[k] = Complex{ spectrum[k].Imaginary(), -spectrum[k].Real() };
+
+        VectorReal& hilbertTd{ fft.Inverse(hilbertSpectrum) };
+
+        for (std::size_t n{ 0 }; n < N; ++n)
+            analytic[n] = Complex{ x[n], hilbertTd[n] };
+
+        return analytic;
+    }
+
+    template<typename T, std::size_t N>
+    T AnalyticSignalFft<T, N>::InstantaneousAmplitude(Complex a)
+    {
+        return std::sqrt(a.Real() * a.Real() + a.Imaginary() * a.Imaginary());
+    }
+
+    template<typename T, std::size_t N>
+    T AnalyticSignalFft<T, N>::InstantaneousPhase(Complex a)
+    {
+        return std::atan2(a.Imaginary(), a.Real());
+    }
+
+    template<typename T, std::size_t N>
+    T AnalyticSignalFft<T, N>::InstantaneousFrequency(T phaseNow, T phasePrev, T ts)
+    {
+        T dphi{ phaseNow - phasePrev };
+        while (dphi > std::numbers::pi_v<T>)
+            dphi -= T(2) * std::numbers::pi_v<T>;
+        while (dphi < -std::numbers::pi_v<T>)
+            dphi += T(2) * std::numbers::pi_v<T>;
+        return dphi / (T(2) * std::numbers::pi_v<T> * ts);
+    }
+
+    /// HilbertFir implementation ///
+
+    template<typename T, std::size_t Taps>
+    HilbertFir<T, Taps>::HilbertFir()
+    {
+        coeff.fill(T(0));
+        delay.fill(T(0));
+        for (std::size_t i{ 0 }; i < Taps; ++i)
+        {
+            int k{ static_cast<int>(i) - static_cast<int>(centerTap) };
+            if (k != 0 && (k % 2) != 0)
+            {
+                T w{ T(0.54) - T(0.46) * std::cos(T(2) * std::numbers::pi_v<T> * static_cast<T>(i) / static_cast<T>(Taps - 1)) };
+                coeff[i] = (T(2) / (std::numbers::pi_v<T> * static_cast<T>(k))) * w;
+            }
+        }
+    }
+
+    template<typename T, std::size_t Taps>
+    OPTIMIZE_FOR_SPEED typename HilbertFir<T, Taps>::Complex HilbertFir<T, Taps>::Filter(T x)
+    {
+        delay[writeIndex] = x;
+
+        T imag{ T(0) };
+        for (std::size_t i{ 0 }; i < Taps; ++i)
+        {
+            std::size_t idx{ (writeIndex + Taps - i) % Taps };
+            imag += coeff[i] * delay[idx];
+        }
+
+        std::size_t realIdx{ (writeIndex + Taps - centerTap) % Taps };
+        T real{ delay[realIdx] };
+
+        writeIndex = (writeIndex + 1) % Taps;
+
+        return Complex{ real, imag };
+    }
+
+#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
+    extern template class AnalyticSignalFft<float, 64>;
+    extern template class HilbertFir<float, 31>;
+#endif
+}

numerical/analysis/test/CMakeLists.txt

@@ -14,6 +14,7 @@ target_sources(numerical.analysis_test PRIVATE
     TestDiscreteCosineTransform.cpp
     TestFastFourierTransformRadix2Impl.cpp
     TestGoertzelAlgorithm.cpp
+    TestHilbertTransform.cpp
     TestPowerDensitySpectrum.cpp
     TestRealFastFourierTransform.cpp
     TestSignalDetectors.cpp