feat(analysis): accept the bin phasor in GoertzelAlgorithm (#360)

A caller can supply e^{jω} of the bin instead of the bin index. The resonator
needs only cos ω and sin ω, evaluated once; on a Cortex-M4 the software
sinf/cosf behind math::Sin/math::Cos pull in several kilobytes of argument
reduction, which a caller with its own trigonometry can now keep out of the
image. The bin-index constructor delegates to the new one and evaluates the
cosine once instead of twice.


Claude-Session: https://claude.ai/code/session_01SeegZiU2Up4ZYQBwnfrLeP

Co-authored-by: Claude <noreply@anthropic.com>

Author: gabrielfrasantos

doc/analysis/GoertzelAlgorithm.md

@@ -90,6 +90,8 @@ The DC bin ($k = 0$) and the Nyquist bin ($k = N/2$) are purely real; their imag
 
 The magnitude-squared form $s_1^2 + s_2^2 - c\,s_1 s_2$ avoids the square root when only power detection (threshold comparison) is needed, saving a transcendental call on platforms without hardware `sqrt`.
 
+The recurrence and the final step need only $\cos\omega$ and $\sin\omega$ of the bin, $\omega = 2\pi k/N$, evaluated once per configuration. A caller that already evaluates them by other means — a sine table refined by a short series, or an exact value for a known angle — can supply the unit bin phasor $e^{j\omega}$ directly. On targets whose sine and cosine come from a software math library, that keeps the library's general argument-reduction code, several kilobytes on a Cortex-M4, out of the image. The phasor's accuracy sets the resonator frequency: an error $\varepsilon$ in $\cos\omega$ detunes it by about $\varepsilon / \sin\omega$ radians per sample, which over a block of $N$ samples rotates the result by about $N\varepsilon / (2\sin\omega)$.
+
 When several tones must be tracked simultaneously, one independent Goertzel instance per tone runs in $O(BN)$ total with $O(B)$ state — still far cheaper than an FFT for small $B$.
 
 A sliding Goertzel variant updates one sample at a time by subtracting the oldest sample, giving a continuous output at the cost of maintaining an input ring buffer; this trades $O(1)$ memory for $O(N)$ memory and is outside the scope of the no-heap constraint here.

numerical/analysis/GoertzelAlgorithm.hpp

@@ -19,6 +19,7 @@ namespace analysis
     public:
         GoertzelAlgorithm(std::size_t k, std::size_t blockLength);
         GoertzelAlgorithm(T targetHz, T sampleHz, std::size_t blockLength);
+        GoertzelAlgorithm(math::Complex<T> binPhasor, std::size_t blockLength);
 
         [[nodiscard]] static std::optional<GoertzelAlgorithm> TryCreate(std::size_t k, std::size_t blockLength);
 
@@ -31,6 +32,8 @@ namespace analysis
         static T Coefficient(std::size_t k, std::size_t blockLength);
 
     private:
+        static math::Complex<T> BinPhasor(std::size_t k, std::size_t blockLength);
+
         T coeff;
         T cosine;
         T sine;
@@ -51,9 +54,14 @@ namespace analysis
 
     template<typename T>
     GoertzelAlgorithm<T>::GoertzelAlgorithm(std::size_t k, std::size_t blockLength)
-        : coeff{ T{ 2 } * math::Cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
-        , cosine{ math::Cos(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
-        , sine{ math::Sin(T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength)) }
+        : GoertzelAlgorithm(BinPhasor(k, blockLength), blockLength)
+    {}
+
+    template<typename T>
+    GoertzelAlgorithm<T>::GoertzelAlgorithm(math::Complex<T> binPhasor, std::size_t blockLength)
+        : coeff{ T{ 2 } * binPhasor.Real() }
+        , cosine{ binPhasor.Real() }
+        , sine{ binPhasor.Imaginary() }
         , blockSize{ blockLength }
     {
         really_assert(blockLength > 0);
@@ -106,6 +114,13 @@ namespace analysis
         sampleCount = 0;
     }
 
+    template<typename T>
+    math::Complex<T> GoertzelAlgorithm<T>::BinPhasor(std::size_t k, std::size_t blockLength)
+    {
+        const T angle{ T{ 2 } * std::numbers::pi_v<T> * static_cast<T>(k) / static_cast<T>(blockLength) };
+        return math::Complex<T>{ math::Cos(angle), math::Sin(angle) };
+    }
+
     template<typename T>
     T GoertzelAlgorithm<T>::Coefficient(std::size_t k, std::size_t blockLength)
     {

numerical/analysis/test/TestGoertzelAlgorithm.cpp

@@ -99,6 +99,32 @@ TEST_F(TestGoertzelAlgorithm, result_matches_direct_dft_bin)
     EXPECT_NEAR(result.Imaginary(), refImag, 1e-3f);
 }
 
+TEST_F(TestGoertzelAlgorithm, supplied_bin_phasor_matches_direct_dft_bin)
+{
+    constexpr std::size_t k{ 3 };
+    const float binAngle{ 2.0f * std::numbers::pi_v<float> * static_cast<float>(k) / static_cast<float>(N) };
+    analysis::GoertzelAlgorithm<float> goertzel{ math::Complex<float>{ std::cos(binAngle), std::sin(binAngle) }, N };
+
+    std::array<float, N> signal{};
+    for (std::size_t n = 0; n < N; ++n)
+        signal[n] = std::sin(binAngle * static_cast<float>(n)) + 0.5f;
+
+    for (float s : signal)
+        goertzel.Push(s);
+
+    float refReal{ 0.0f };
+    float refImag{ 0.0f };
+    for (std::size_t n = 0; n < N; ++n)
+    {
+        refReal += signal[n] * std::cos(binAngle * static_cast<float>(n));
+        refImag -= signal[n] * std::sin(binAngle * static_cast<float>(n));
+    }
+
+    math::Complex<float> result{ goertzel.Result() };
+    EXPECT_NEAR(result.Real(), refReal, 1e-3f);
+    EXPECT_NEAR(result.Imaginary(), refImag, 1e-3f);
+}
+
 TEST_F(TestGoertzelAlgorithm, magnitude_squared_equals_result_norm_squared)
 {
     constexpr std::size_t k{ 2 };