Merge branch 'main' into feature/add-goertzel-algo

Author: gabrielfrasantos

ROADMAP.md

@@ -28,7 +28,6 @@ Difficulty legend:
 | #  | Component                                            | Target module             | Difficulty |
 |----|------------------------------------------------------|---------------------------|------------|
 | 10 | Gain-scheduled controller                            | `controllers`             | ★★☆☆☆      |
-| 11 | Convolution & correlation utilities                  | `analysis`                | ★★☆☆☆      |
 | 12 | Polynomial least-squares curve fitting               | `estimators/offline`      | ★★☆☆☆      |
 | 14 | CIC (Cascaded Integrator-Comb) filter                | `filters/passive`         | ★★☆☆☆      |
 | 15 | Biquad / Second-Order-Section cascade                | `filters/passive`         | ★★★☆☆      |

doc/analysis/ConvolutionCorrelation.md

@@ -0,0 +1,145 @@
+# Convolution & Correlation
+
+## Overview & Motivation
+
+Two signals interacting in time — one being filtered, delayed, or matched against another — are
+described by convolution and correlation. Convolution is the time-domain action of any linear
+time-invariant (LTI) filter: given an input signal and a system's impulse response, their
+convolution is the output. Correlation measures the similarity between two signals as a function
+of relative shift (lag), making it the foundation of matched filtering, pulse detection, and
+time-delay estimation.
+
+Both operations are indispensable in embedded DSP: they underlie FIR filtering, pitch detection,
+radar/sonar processing, and auto-regressive modeling. Bounded, no-heap implementations make
+these algorithms practical on microcontrollers.
+
+## Mathematical Theory
+
+### Linear Convolution
+
+For finite sequences $x$ of length $M$ and $h$ of length $K$, the linear convolution $y = x * h$
+is a sequence of length $M + K - 1$ defined by:
+
+$$y[n] = \sum_{k=0}^{M-1} x[k] \cdot h[n - k], \quad n = 0, \ldots, M+K-2$$
+
+where $h[j] = 0$ for $j < 0$ or $j \geq K$.
+
+### Circular Convolution
+
+For two sequences of equal length $N$, circular (periodic) convolution wraps indices modulo $N$:
+
+$$y[n] = \sum_{k=0}^{N-1} x[k] \cdot h[(n - k) \bmod N], \quad n = 0, \ldots, N-1$$
+
+Circular convolution equals linear convolution only when both operands are zero-padded to at least
+$M + K - 1$ samples. Without padding, time-aliasing occurs.
+
+### Cross-Correlation
+
+Cross-correlation of $x$ (length $M$) and $y$ (length $K$) is defined as:
+
+$$r_{xy}[\ell] = \sum_{n} x[n] \cdot y[n - \ell]$$
+
+This equals linear convolution of $x$ with the time-reversed $y$:
+
+$$r_{xy} = x * \overline{y}$$
+
+### Auto-Correlation
+
+Auto-correlation is cross-correlation of a signal with itself:
+
+$$r_{xx}[\ell] = \sum_{n} x[n] \cdot x[n - \ell]$$
+
+It is always symmetric about the zero-lag index, and its maximum occurs at zero lag.
+Dividing by $r_{xx}[0]$ yields a normalized coefficient in $[-1, 1]$.
+
+### Convolution Theorem (Fast Path)
+
+For long signals, the convolution theorem enables an $O(L \log L)$ computation:
+
+$$x * h = \mathcal{F}^{-1}\!\left(\mathcal{F}(x) \cdot \mathcal{F}(h)\right)$$
+
+where $\mathcal{F}$ denotes the DFT and $L$ is the next power of two $\geq M + K - 1$. The
+pointwise complex product replaces the $O(MK)$ direct sum.
+
+## Complexity Analysis
+
+| Operation            | Time          | Space (extra) | Notes                                 |
+|----------------------|---------------|---------------|---------------------------------------|
+| Linear convolution   | $O(MK)$       | $O(1)$        | Writes into caller-owned buffer       |
+| Circular convolution | $O(N^2)$      | $O(1)$        | Same as linear for equal-length input |
+| Cross-correlation    | $O(MK)$       | $O(K)$        | Reverses one operand on the stack     |
+| Auto-correlation     | $O(M^2)$      | $O(1)$        | Alias of cross-correlation            |
+| Fast convolution     | $O(L \log L)$ | $O(L)$        | $L = 2^{\lceil\log_2(M+K-1)\rceil}$   |
+
+The fast path is beneficial when $MK > L \log_2 L$, roughly when both operands exceed 32–64
+samples.
+
+## Step-by-Step Walkthrough
+
+Linear convolution of $x = [1, 2, 3]$ and $h = [1, 1]$ (output length $= 4$):
+
+| $n$ | Active $k$ range | Computation             | $y[n]$ |
+|-----|------------------|-------------------------|--------|
+| 0   | $k=0$            | $1 \cdot 1$             | 1      |
+| 1   | $k=0,1$          | $2 \cdot 1 + 1 \cdot 1$ | 3      |
+| 2   | $k=1,2$          | $3 \cdot 1 + 2 \cdot 1$ | 5      |
+| 3   | $k=2$            | $3 \cdot 1$             | 3      |
+
+Result: $[1, 3, 5, 3]$.
+
+Auto-correlation of $[1, 1, 1, 1]$ produces the triangular sequence $[1, 2, 3, 4, 3, 2, 1]$,
+with the peak at the zero-lag center index (index 3 in 0-based notation).
+
+## Pitfalls & Edge Cases
+
+- **Time-aliasing**: Circular convolution of length $N$ equals linear convolution only when both
+  operands are zero-padded to $M + K - 1 \leq N$. Failing to pad causes energy from the tail to
+  wrap around and corrupt the beginning of the output.
+- **Lag estimation**: `ArgMaxLag` returns the index of the maximum in the cross-correlation
+  output. For cross-correlation of two length-$M$ sequences the zero-lag corresponds to index
+  $M - 1$; a delay of $d$ samples appears at index $M - 1 + d$.
+- **Normalization**: The raw auto-correlation at zero lag equals the signal energy. Dividing the
+  entire sequence by $r[0]$ gives a coefficient bounded in $[-1, 1]$.
+- **Float precision**: Accumulating many products in single precision can lose significance for
+  very long sequences. The MAC loop uses the native float accumulator; for sequences longer than
+  a few hundred samples, extended precision may be warranted.
+
+## Variants & Generalizations
+
+- **Normalized cross-correlation**: Divide $r_{xy}[\ell]$ by $\sqrt{r_{xx}[0] \cdot r_{yy}[0]}$
+  for an amplitude-invariant similarity measure.
+- **Overlap-add / overlap-save**: Partition long signals into short blocks and convolve each
+  with $h$ via the FFT fast path, then reassemble. Reduces latency for streaming applications.
+- **Partial correlation**: Compute cross-correlation over a sliding window rather than the full
+  signal length — useful for non-stationary signals.
+
+## Applications
+
+- **FIR filtering**: Direct-form FIR is linear convolution of the input with the filter's impulse
+  response.
+- **Matched filtering**: Cross-correlate the received signal with a known template to detect it
+  in noise.
+- **Time-delay estimation (TDOA)**: The lag of the cross-correlation peak between two sensors
+  gives the inter-sensor arrival time difference.
+- **Pitch detection**: Auto-correlation of a speech frame reveals the fundamental period as the
+  first large secondary peak.
+- **AR modeling**: Yule-Walker equations use auto-correlation values to fit auto-regressive
+  parameters.
+
+## Connections to Other Algorithms
+
+- **FastFourierTransform**: The fast convolution path is a direct application of the FFT-based
+  convolution theorem.
+- **GoertzelAlgorithm**: A single-frequency DFT bin — cross-correlating with a complex
+  exponential at one frequency is equivalent.
+- **YuleWalker**: Constructs the AR system matrix from auto-correlation lags computed by
+  `AutoCorrelation`.
+- **PowerSpectralDensity**: The power spectrum is the DFT of the auto-correlation sequence
+  (Wiener–Khinchin theorem).
+
+## References & Further Reading
+
+- A. V. Oppenheim, R. W. Schafer, *Discrete-Time Signal Processing*, 3rd ed., Prentice Hall,
+  2009, Chapters 2 & 8.
+- A. V. Oppenheim, A. S. Willsky, *Signals and Systems*, 2nd ed., Prentice Hall, 1997,
+  Chapter 9.

numerical/analysis/CMakeLists.txt

@@ -12,6 +12,7 @@ target_link_libraries(numerical.analysis ${NUMERICAL_VISIBILITY}
 )
 
 target_sources(numerical.analysis PRIVATE
+    ConvolutionCorrelation.hpp
     DiscreteCosineTransform.hpp
     FastFourierTransform.hpp
     FastFourierTransformRadix2Impl.hpp
@@ -21,6 +22,7 @@ target_sources(numerical.analysis PRIVATE
 )
 
 numerical_add_coverage_sources(numerical.analysis
+    ConvolutionCorrelation.cpp
     FastFourierTransformRadix2Impl.cpp
     GoertzelAlgorithm.cpp
     PowerDensitySpectrum.cpp

numerical/analysis/ConvolutionCorrelation.cpp

@@ -0,0 +1,32 @@
+#include "numerical/analysis/ConvolutionCorrelation.hpp"
+
+namespace analysis
+{
+    template void LinearConvolution<float, 3, 3>(
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        infra::BoundedVector<float>::WithMaxSize<5>&);
+
+    template void CircularConvolution<float, 4>(
+        const infra::BoundedVector<float>::WithMaxSize<4>&,
+        const infra::BoundedVector<float>::WithMaxSize<4>&,
+        infra::BoundedVector<float>::WithMaxSize<4>&);
+
+    template void CrossCorrelation<float, 5, 5>(
+        const infra::BoundedVector<float>::WithMaxSize<5>&,
+        const infra::BoundedVector<float>::WithMaxSize<5>&,
+        infra::BoundedVector<float>::WithMaxSize<9>&);
+
+    template void AutoCorrelation<float, 4>(
+        const infra::BoundedVector<float>::WithMaxSize<4>&,
+        infra::BoundedVector<float>::WithMaxSize<7>&);
+
+    template std::size_t ArgMaxLag<float, 9>(
+        const infra::BoundedVector<float>::WithMaxSize<9>&);
+
+    template void FastConvolution<float, 3, 3, 8>(
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        infra::BoundedVector<float>::WithMaxSize<5>&,
+        FastFourierTransform<float>&);
+}

numerical/analysis/ConvolutionCorrelation.hpp

@@ -0,0 +1,163 @@
+#pragma once
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC optimize("O3", "fast-math")
+#endif
+
+#include "numerical/analysis/FastFourierTransform.hpp"
+#include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/ComplexNumber.hpp"
+#include <cstddef>
+#include <infra/util/BoundedVector.hpp>
+#include <type_traits>
+
+namespace analysis
+{
+    template<typename T, std::size_t M, std::size_t K>
+    OPTIMIZE_FOR_SPEED void LinearConvolution(
+        const typename infra::BoundedVector<T>::template WithMaxSize<M>& x,
+        const typename infra::BoundedVector<T>::template WithMaxSize<K>& h,
+        typename infra::BoundedVector<T>::template WithMaxSize<M + K - 1>& y)
+    {
+        static_assert(std::is_floating_point_v<T>, "ConvolutionCorrelation supports floating-point types only");
+        y.clear();
+        y.resize(M + K - 1, T{ 0 });
+        for (std::size_t n = 0; n < M + K - 1; ++n)
+        {
+            T acc{ 0 };
+            std::size_t kStart{ n >= K - 1 ? n - (K - 1) : 0 };
+            std::size_t kEnd{ n < M - 1 ? n : M - 1 };
+            for (std::size_t k = kStart; k <= kEnd; ++k)
+                acc += x[k] * h[n - k];
+            y[n] = acc;
+        }
+    }
+
+    template<typename T, std::size_t N>
+    OPTIMIZE_FOR_SPEED void CircularConvolution(
+        const typename infra::BoundedVector<T>::template WithMaxSize<N>& x,
+        const typename infra::BoundedVector<T>::template WithMaxSize<N>& h,
+        typename infra::BoundedVector<T>::template WithMaxSize<N>& y)
+    {
+        static_assert(std::is_floating_point_v<T>, "ConvolutionCorrelation supports floating-point types only");
+        y.clear();
+        y.resize(N, T{ 0 });
+        for (std::size_t n = 0; n < N; ++n)
+        {
+            T acc{ 0 };
+            for (std::size_t k = 0; k < N; ++k)
+                acc += x[k] * h[(n + N - k) % N];
+            y[n] = acc;
+        }
+    }
+
+    template<typename T, std::size_t M, std::size_t K>
+    OPTIMIZE_FOR_SPEED void CrossCorrelation(
+        const typename infra::BoundedVector<T>::template WithMaxSize<M>& x,
+        const typename infra::BoundedVector<T>::template WithMaxSize<K>& y,
+        typename infra::BoundedVector<T>::template WithMaxSize<M + K - 1>& r)
+    {
+        static_assert(std::is_floating_point_v<T>, "ConvolutionCorrelation supports floating-point types only");
+        typename infra::BoundedVector<T>::template WithMaxSize<K> yRev;
+        yRev.resize(K, T{ 0 });
+        for (std::size_t i = 0; i < K; ++i)
+            yRev[i] = y[K - 1 - i];
+        LinearConvolution<T, M, K>(x, yRev, r);
+    }
+
+    template<typename T, std::size_t M>
+    void AutoCorrelation(
+        const typename infra::BoundedVector<T>::template WithMaxSize<M>& x,
+        typename infra::BoundedVector<T>::template WithMaxSize<2 * M - 1>& r)
+    {
+        static_assert(std::is_floating_point_v<T>, "ConvolutionCorrelation supports floating-point types only");
+        CrossCorrelation<T, M, M>(x, x, r);
+    }
+
+    template<typename T, std::size_t N>
+    std::size_t ArgMaxLag(const typename infra::BoundedVector<T>::template WithMaxSize<N>& r)
+    {
+        static_assert(std::is_floating_point_v<T>, "ConvolutionCorrelation supports floating-point types only");
+        std::size_t best{ 0 };
+        for (std::size_t i = 1; i < N; ++i)
+            if (r[i] > r[best])
+                best = i;
+        return best;
+    }
+
+    template<typename T, std::size_t M, std::size_t K, std::size_t L>
+    void FastConvolution(
+        const typename infra::BoundedVector<T>::template WithMaxSize<M>& x,
+        const typename infra::BoundedVector<T>::template WithMaxSize<K>& h,
+        typename infra::BoundedVector<T>::template WithMaxSize<M + K - 1>& y,
+        FastFourierTransform<T>& fft)
+    {
+        static_assert(std::is_floating_point_v<T>, "ConvolutionCorrelation supports floating-point types only");
+        static_assert((L & (L - 1)) == 0, "L must be a power of two");
+        static_assert(L >= M + K - 1, "L must be at least M+K-1");
+
+        typename infra::BoundedVector<T>::template WithMaxSize<L> xPad;
+        typename infra::BoundedVector<T>::template WithMaxSize<L> hPad;
+        xPad.resize(L, T{ 0 });
+        hPad.resize(L, T{ 0 });
+
+        for (std::size_t i = 0; i < M; ++i)
+            xPad[i] = x[i];
+        for (std::size_t i = 0; i < K; ++i)
+            hPad[i] = h[i];
+
+        typename infra::BoundedVector<math::Complex<T>>::template WithMaxSize<L> X;
+        X.resize(L, math::Complex<T>{ T{ 0 }, T{ 0 } });
+        {
+            auto& rawX = fft.Forward(xPad);
+            for (std::size_t i = 0; i < L; ++i)
+                X[i] = rawX[i];
+        }
+        auto& H = fft.Forward(hPad);
+
+        typename infra::BoundedVector<math::Complex<T>>::template WithMaxSize<L> Y;
+        Y.resize(L, math::Complex<T>{ T{ 0 }, T{ 0 } });
+        for (std::size_t i = 0; i < L; ++i)
+            Y[i] = X[i] * H[i];
+
+        auto& yFull = fft.Inverse(Y);
+
+        y.clear();
+        y.resize(M + K - 1, T{ 0 });
+        for (std::size_t i = 0; i < M + K - 1; ++i)
+            y[i] = yFull[i];
+    }
+}
+
+#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
+namespace analysis
+{
+    extern template void LinearConvolution<float, 3, 3>(
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        infra::BoundedVector<float>::WithMaxSize<5>&);
+
+    extern template void CircularConvolution<float, 4>(
+        const infra::BoundedVector<float>::WithMaxSize<4>&,
+        const infra::BoundedVector<float>::WithMaxSize<4>&,
+        infra::BoundedVector<float>::WithMaxSize<4>&);
+
+    extern template void CrossCorrelation<float, 5, 5>(
+        const infra::BoundedVector<float>::WithMaxSize<5>&,
+        const infra::BoundedVector<float>::WithMaxSize<5>&,
+        infra::BoundedVector<float>::WithMaxSize<9>&);
+
+    extern template void AutoCorrelation<float, 4>(
+        const infra::BoundedVector<float>::WithMaxSize<4>&,
+        infra::BoundedVector<float>::WithMaxSize<7>&);
+
+    extern template std::size_t ArgMaxLag<float, 9>(
+        const infra::BoundedVector<float>::WithMaxSize<9>&);
+
+    extern template void FastConvolution<float, 3, 3, 8>(
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        const infra::BoundedVector<float>::WithMaxSize<3>&,
+        infra::BoundedVector<float>::WithMaxSize<5>&,
+        FastFourierTransform<float>&);
+}
+#endif

numerical/analysis/test/CMakeLists.txt

@@ -9,6 +9,7 @@ target_link_libraries(numerical.analysis_test PUBLIC
 )
 
 target_sources(numerical.analysis_test PRIVATE
+    TestConvolutionCorrelation.cpp
     TestDiscreteCosineTransform.cpp
     TestFastFourierTransformRadix2Impl.cpp
     TestGoertzelAlgorithm.cpp