feat: improving analysis subset

doc/analysis/GoertzelAlgorithm.md

@@ -22,15 +22,15 @@ The DFT sum can be reorganised as the output of a single-pole IIR filter evaluat
 
 $$s[n] = x[n] + c\,s[n-1] - s[n-2], \quad s[-1] = s[-2] = 0$$
 
-After $N$ steps the DFT bin follows from the final two state values:
+After the $N$ input samples, one further step with $x[N] = 0$ gives $s[N] = c\,s[N-1] - s[N-2]$, and the DFT bin follows:
 
-$$X[k] = s[N-1] - s[N-2]\,e^{-j2\pi k/N}$$
+$$X[k] = s[N] - s[N-1]\,e^{-j2\pi k/N}$$
 
 Separating real and imaginary parts:
 
-$$\mathrm{Re}\{X[k]\} = s[N-1] - s[N-2]\cos\!\left(\tfrac{2\pi k}{N}\right)$$
+$$\mathrm{Re}\{X[k]\} = s[N] - s[N-1]\cos\!\left(\tfrac{2\pi k}{N}\right)$$
 
-$$\mathrm{Im}\{X[k]\} = s[N-2]\sin\!\left(\tfrac{2\pi k}{N}\right)$$
+$$\mathrm{Im}\{X[k]\} = s[N-1]\sin\!\left(\tfrac{2\pi k}{N}\right)$$
 
 The magnitude can be obtained without the final trigonometric products using the identity:
 
@@ -64,17 +64,19 @@ $c = 2\cos(\pi/2) = 0$
 | $n$ | $x[n]$ | $s[n] = x[n] + 0\cdot s[n-1] - s[n-2]$ |
 |-----|--------|----------------------------------------|
 | 0   | 1      | $1 + 0 - 0 = 1$                        |
-| 1   | 0      | $0 + 0 - 1 = -1$                       |
-| 2   | -1     | $-1 + 0 - 0 = -1$  (note: $s[-1]=0$)   |
-| 3   | 0      | $0 + 0 - (-1) = 1$                     |
+| 1   | 0      | $0 + 0 - 0 = 0$  (note: $s[-1]=0$)     |
+| 2   | -1     | $-1 + 0 - 1 = -2$                      |
+| 3   | 0      | $0 + 0 - 0 = 0$                        |
+
+Extra step: $s[4] = 0\cdot s[3] - s[2] = -(-2) = 2$.
 
 $\cos(2\pi/4) = 0$, $\sin(2\pi/4) = 1$
 
-$\mathrm{Re}\{X[1]\} = 1 - (-1)\cdot 0 = 1$
+$\mathrm{Re}\{X[1]\} = s[4] - s[3]\cdot 0 = 2$
 
-$\mathrm{Im}\{X[1]\} = (-1)\cdot 1 = -1$
+$\mathrm{Im}\{X[1]\} = s[3]\cdot 1 = 0$
 
-Direct DFT check: $X[1] = 1 - j$.  Matches.
+Direct DFT check: $X[1] = 1 - (-1) = 2$.  Matches.
 
 ## Pitfalls & Edge Cases
 

numerical/analysis/DiscreteCosineTransform.hpp

@@ -31,6 +31,7 @@ namespace analysis
     private:
         FastFourierTransform<QNumberType>& fft;
         typename infra::BoundedVector<QNumberType>::template WithMaxSize<Length> output;
+        typename infra::BoundedVector<QNumberType>::template WithMaxSize<Length> reordered;
         typename VectorComplex::template WithMaxSize<Length> complexBuffer;
     };
 
@@ -41,6 +42,7 @@ namespace analysis
         : fft(fft)
     {
         output.resize(Length);
+        reordered.resize(Length);
         complexBuffer.resize(Length);
     }
 
@@ -49,7 +51,13 @@ namespace analysis
         typename DiscreteConsineTransform<QNumberType, Length>::VectorReal&
         DiscreteConsineTransform<QNumberType, Length>::Forward(VectorReal& input)
     {
-        auto& fftResult = fft.Forward(input);
+        for (std::size_t n = 0; n < Length / 2; ++n)
+        {
+            reordered[n] = input[2 * n];
+            reordered[Length - 1 - n] = input[2 * n + 1];
+        }
+
+        auto& fftResult = fft.Forward(reordered);
 
         output[0] = QNumberType(math::ToFloat(fftResult[0].Real()) / std::sqrt(static_cast<float>(Length)));
 
@@ -70,18 +78,31 @@ namespace analysis
     template<typename QNumberType, std::size_t Length>
     typename DiscreteConsineTransform<QNumberType, Length>::VectorReal& DiscreteConsineTransform<QNumberType, Length>::Inverse(VectorReal& input)
     {
-        complexBuffer[0] = math::Complex<QNumberType>{ QNumberType(math::ToFloat(input[0]) * std::sqrt(static_cast<float>(Length))), QNumberType(0.0f) };
+        float sqrtN = std::sqrt(static_cast<float>(Length));
+
+        complexBuffer[0] = math::Complex<QNumberType>{ QNumberType(math::ToFloat(input[0]) * sqrtN), QNumberType(0.0f) };
 
         for (std::size_t k = 1; k < Length; ++k)
         {
+            float real = math::ToFloat(input[k]) * sqrtN / 2.0f;
+            float imag = -math::ToFloat(input[Length - k]) * sqrtN / 2.0f;
+
             float angle = static_cast<float>(k) * std::numbers::pi_v<float> / (2.0f * static_cast<float>(Length));
-            float scale = std::sqrt(static_cast<float>(Length)) / 2.0f;
+            float cosine = std::cos(angle);
+            float sine = std::sin(angle);
 
-            float value = math::ToFloat(input[k]) * scale;
-            complexBuffer[k] = math::Complex<QNumberType>{ QNumberType(value * std::cos(angle)), QNumberType(value * std::sin(angle)) };
+            complexBuffer[k] = math::Complex<QNumberType>{ QNumberType(real * cosine - imag * sine), QNumberType(real * sine + imag * cosine) };
         }
 
-        return fft.Inverse(complexBuffer);
+        auto& timeDomain = fft.Inverse(complexBuffer);
+
+        for (std::size_t n = 0; n < Length / 2; ++n)
+        {
+            output[2 * n] = timeDomain[n];
+            output[2 * n + 1] = timeDomain[Length - 1 - n];
+        }
+
+        return output;
     }
 
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD

numerical/analysis/GoertzelAlgorithm.hpp

@@ -73,8 +73,9 @@ namespace analysis
     template<typename T>
     math::Complex<T> GoertzelAlgorithm<T>::Result() const
     {
-        T real{ s1 - s2 * cosine };
-        T imag{ s2 * sine };
+        T sExtra{ coeff * s1 - s2 };
+        T real{ sExtra - s1 * cosine };
+        T imag{ s1 * sine };
         return math::Complex<T>{ real, imag };
     }
 

numerical/analysis/PowerDensitySpectrum.cpp

@@ -4,4 +4,5 @@
 namespace analysis
 {
     template class PowerSpectralDensity<float, 512, test::FftStub<float, 512>, test::TwiddleFactorsStub<float, 256>, 50>;
+    template class PowerSpectralDensity<float, 512, test::FftStub<float, 512>, test::TwiddleFactorsStub<float, 256>, 0>;
 }

numerical/analysis/PowerDensitySpectrum.hpp

@@ -46,10 +46,8 @@ namespace analysis
         void ResetSegment();
 
     private:
-        static constexpr QNumberType segmentSizeInverted = QNumberType(1.0f / static_cast<float>(SegmentSize));
         static constexpr std::size_t overlapSize = (SegmentSize * Overlap) / 100;
         static constexpr std::size_t step = SegmentSize - overlapSize;
-        QNumberType frequencyResolution;
         windowing::Window<QNumberType>& window;
         QNumberType samplingTimeInSeconds;
         TwiddleFactor twiddleFactors;
@@ -65,7 +63,6 @@ namespace analysis
         windowing::Window<QNumberType>& window, QNumberType samplingTimeInSeconds)
         : window(window)
         , samplingTimeInSeconds(samplingTimeInSeconds)
-        , frequencyResolution(QNumberType(samplingTimeInSeconds * segmentSizeInverted))
     {}
 
     template<typename QNumberType, std::size_t SegmentSize, typename Fft, typename TwiddleFactor, std::size_t Overlap>
@@ -76,6 +73,8 @@ namespace analysis
 
         ResetOutput();
 
+        std::size_t segmentCount = 0;
+
         for (std::size_t i = 0; i + SegmentSize <= input.size(); i += step)
         {
             ResetSegment();
@@ -87,10 +86,15 @@ namespace analysis
 
             for (std::size_t k = 0; k <= SegmentSize / 2; ++k)
                 y[k] += QNumberType(math::ToFloat(MagnitudeSquared(spectrum[k])) / static_cast<float>(SegmentSize));
+
+            ++segmentCount;
         }
 
+        float normalization = math::ToFloat(samplingTimeInSeconds) /
+                              (math::ToFloat(window.Power(SegmentSize)) * static_cast<float>(segmentCount));
+
         for (std::size_t i = 0; i < y.size(); ++i)
-            y[i] *= frequencyResolution;
+            y[i] = QNumberType(math::ToFloat(y[i]) * normalization);
 
         return y;
     }
@@ -118,5 +122,6 @@ namespace analysis
 
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
     extern template class PowerSpectralDensity<float, 512, test::FftStub<float, 512>, test::TwiddleFactorsStub<float, 256>, 50>;
+    extern template class PowerSpectralDensity<float, 512, test::FftStub<float, 512>, test::TwiddleFactorsStub<float, 256>, 0>;
 #endif
 }

numerical/analysis/test/TestDecibels.cpp

@@ -39,6 +39,11 @@ TEST_F(TestDecibels, negative_ratio_returns_floor)
     EXPECT_NEAR(analysis::ToDecibels(-1.0f), analysis::DecibelFloor<float>::value, math::Tolerance<float>());
 }
 
+TEST_F(TestDecibels, tiny_positive_ratio_clamps_to_floor)
+{
+    EXPECT_NEAR(analysis::ToDecibels(1e-9f), analysis::DecibelFloor<float>::value, math::Tolerance<float>());
+}
+
 TEST_F(TestDecibels, from_decibels_twenty_returns_ten)
 {
     EXPECT_NEAR(analysis::FromDecibels(20.0f), 10.0f, math::Tolerance<float>());
@@ -56,5 +61,5 @@ TEST_F(TestDecibels, attenuation_db_computes_difference)
 
 TEST_F(TestDecibels, ripple_db_computes_passband_variation)
 {
-    EXPECT_NEAR(analysis::RippleDb(1.0f, 0.9f), analysis::ToDecibels(1.0f) - analysis::ToDecibels(0.9f), math::Tolerance<float>());
+    EXPECT_NEAR(analysis::RippleDb(1.0f, 0.9f), 0.9151f, 1e-3f);
 }