fix: resolve issues 251, 252, and 253

doc/analysis/DiscreteCosineTransform.md

@@ -10,21 +10,25 @@ This library implements **DCT-II** (the most common variant) by leveraging the F
 
 ### DCT-II Definition
 
-For a length-$N$ sequence $x[n]$, the DCT-II is defined as:
+This implementation uses the **orthonormal** DCT-II. For a length-$N$ sequence $x[n]$:
 
-$$X[k] = \sum_{n=0}^{N-1} x[n] \cos\!\left(\frac{\pi}{N}\left(n + \tfrac{1}{2}\right) k\right), \quad k = 0, 1, \ldots, N-1$$
+$$X[0] = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x[n]$$
 
-The inverse (DCT-III, also called IDCT) recovers $x[n]$:
+$$X[k] = \sqrt{\frac{2}{N}} \sum_{n=0}^{N-1} x[n] \cos\!\left(\frac{\pi}{N}\left(n + \tfrac{1}{2}\right) k\right), \quad k = 1, \ldots, N-1$$
 
-$$x[n] = \frac{1}{N}\left[\frac{X[0]}{2} + \sum_{k=1}^{N-1} X[k] \cos\!\left(\frac{\pi}{N}\left(n + \tfrac{1}{2}\right) k\right)\right]$$
+This normalisation makes the transform matrix unitary, so Parseval's theorem holds exactly: $\sum_k X[k]^2 = \sum_n x[n]^2$. The inverse (orthonormal DCT-III) recovers $x[n]$:
+
+$$x[n] = \frac{X[0]}{\sqrt{N}} + \sqrt{\frac{2}{N}} \sum_{k=1}^{N-1} X[k] \cos\!\left(\frac{\pi}{N}\left(n + \tfrac{1}{2}\right) k\right)$$
 
 ### FFT-Based Computation
 
 The DCT can be computed via the DFT of a reordered sequence:
 
 1. Form a new sequence $y[n]$ by interleaving even-indexed and reverse odd-indexed samples of $x$.
 2. Compute the $N$-point FFT: $Y[k] = \text{FFT}\{y\}$.
-3. Apply twiddle factors: $X[k] = 2 \cdot \text{Re}\!\left(W[k] \cdot Y[k]\right)$, where $W[k] = e^{-j\pi k / 2N}$.
+3. Apply twiddle factors and orthonormal scale:
+   - $k = 0$: $X[0] = \text{Re}(Y[0]) / \sqrt{N}$
+   - $k \ge 1$: $X[k] = \sqrt{2/N} \cdot \text{Re}\!\left(W[k] \cdot Y[k]\right)$, where $W[k] = e^{-j\pi k / 2N}$.
 
 ### Why Cosine Basis?
 
@@ -50,25 +54,25 @@ $$y = [x[0],\; x[2],\; x[3],\; x[1]] = [1, 3, 4, 2]$$
 
 **Step 2 — Compute FFT**
 
-$$Y = \text{FFT}([1, 3, 4, 2]) = [10,\; -3+j,\; -2,\; -3-j]$$
+$$Y = \text{FFT}([1, 3, 4, 2]) = [10,\; -3-j,\; 0,\; -3+j]$$
 
-**Step 3 — Apply twiddle factors** $W[k] = e^{-j\pi k/8}$
+**Step 3 — Apply twiddle factors and orthonormal scale** ($N=4$, $W[k] = e^{-j\pi k/8}$)
 
-| $k$ | $W[k]$         | $W[k] \cdot Y[k]$ | $X[k] = 2 \cdot \text{Re}(\cdot)$ |
-|-----|----------------|-------------------|-----------------------------------|
-| 0   | 1              | 10                | 20                                |
-| 1   | $e^{-j\pi/8}$  | ≈ −2.22 − 1.90j   | ≈ −4.44                           |
-| 2   | $e^{-j\pi/4}$  | ≈ −1.41 + 1.41j   | ≈ −2.83                           |
-| 3   | $e^{-j3\pi/8}$ | ≈ −0.24 + 3.07j   | ≈ −0.47                           |
+| $k$ | $W[k]$         | $W[k] \cdot Y[k]$         | $X[k]$                                    |
+|-----|----------------|---------------------------|-------------------------------------------|
+| 0   | —              | —                         | $10 / \sqrt{4} = 5$                       |
+| 1   | $e^{-j\pi/8}$  | $\approx -3.154 + 0.224j$ | $\sqrt{0.5} \cdot (-3.154) \approx -2.23$ |
+| 2   | $e^{-j\pi/4}$  | $0$                       | $0$                                       |
+| 3   | $e^{-j3\pi/8}$ | $\approx -0.224 + 3.154j$ | $\sqrt{0.5} \cdot (-0.224) \approx -0.16$ |
 
-**Output:** $X \approx [20, -4.44, -2.83, -0.47]$
+**Output:** $X \approx [5,\; -2.23,\; 0,\; -0.16]$
 
-Notice how most of the energy is in $X[0]$ (the DC component) — energy compaction in action.
+Notice how most of the energy is in $X[0]$ (the DC component) — energy compaction in action. The orthonormal scale ensures $\sum_k X[k]^2 = 1^2+2^2+3^2+4^2 = 30$.
 
 ## Pitfalls & Edge Cases
 
 - **Power-of-2 length required** — inherited from the underlying FFT constraint.
-- **Normalization convention.** Different references use different scaling (some include $\sqrt{2/N}$). Verify which convention the consumer expects.
+- **Normalization convention.** This library uses the **orthonormal** convention ($1/\sqrt{N}$ for $k=0$, $\sqrt{2/N}$ for $k \ge 1$), which satisfies Parseval's theorem and makes the transform matrix unitary. Other references may use the unnormalized form; scale accordingly when interfacing.
 - **Fixed-point overflow.** The reordering and FFT steps must preserve range; apply the 0.9999 scaling factor used throughout this library.
 - **Inverse accuracy.** Rounding errors accumulate in the forward-then-inverse round-trip, especially for Q15 types.
 - **Real input only.** Complex inputs are not supported by the reordering trick.

numerical/analysis/DiscreteCosineTransform.cpp

@@ -3,7 +3,7 @@
 
 namespace analysis
 {
-    template class DiscreteConsineTransform<float, 8>;
-    template class DiscreteConsineTransform<math::Q15, 8>;
-    template class DiscreteConsineTransform<math::Q31, 8>;
+    template class DiscreteCosineTransform<float, 8>;
+    template class DiscreteCosineTransform<math::Q15, 8>;
+    template class DiscreteCosineTransform<math::Q31, 8>;
 }

numerical/analysis/DiscreteCosineTransform.hpp

@@ -14,17 +14,17 @@
 namespace analysis
 {
     template<typename QNumberType, std::size_t Length>
-    class DiscreteConsineTransform
+    class DiscreteCosineTransform
     {
-        static_assert((Length & (Length - 1)) == 0, "DiscreteConsineTransform size must be a power of 2");
+        static_assert((Length & (Length - 1)) == 0, "DiscreteCosineTransform size must be a power of 2");
         static_assert(math::is_qnumber<QNumberType>::value || std::is_floating_point_v<QNumberType>,
-            "DiscreteConsineTransform can only be instantiated with math::QNumber types or floating point.");
+            "DiscreteCosineTransform can only be instantiated with math::QNumber types or floating point.");
 
     public:
         using VectorReal = typename FastFourierTransform<QNumberType>::VectorReal;
         using VectorComplex = typename FastFourierTransform<QNumberType>::VectorComplex;
 
-        explicit DiscreteConsineTransform(FastFourierTransform<QNumberType>& fft);
+        explicit DiscreteCosineTransform(FastFourierTransform<QNumberType>& fft);
         VectorReal& Forward(VectorReal& input);
         VectorReal& Inverse(VectorReal& input);
 
@@ -38,7 +38,7 @@ namespace analysis
     // Implementation //
 
     template<typename QNumberType, std::size_t Length>
-    DiscreteConsineTransform<QNumberType, Length>::DiscreteConsineTransform(FastFourierTransform<QNumberType>& fft)
+    DiscreteCosineTransform<QNumberType, Length>::DiscreteCosineTransform(FastFourierTransform<QNumberType>& fft)
         : fft(fft)
     {
         output.resize(Length);
@@ -48,8 +48,8 @@ namespace analysis
 
     template<typename QNumberType, std::size_t Length>
     OPTIMIZE_FOR_SPEED
-        typename DiscreteConsineTransform<QNumberType, Length>::VectorReal&
-        DiscreteConsineTransform<QNumberType, Length>::Forward(VectorReal& input)
+        typename DiscreteCosineTransform<QNumberType, Length>::VectorReal&
+        DiscreteCosineTransform<QNumberType, Length>::Forward(VectorReal& input)
     {
         for (std::size_t n = 0; n < Length / 2; ++n)
         {
@@ -64,7 +64,7 @@ namespace analysis
         for (std::size_t k = 1; k < Length; ++k)
         {
             float angle = -static_cast<float>(k) * std::numbers::pi_v<float> / (2.0f * Length);
-            float scale = 2.0f / math::Sqrt(static_cast<float>(Length));
+            float scale = math::Sqrt(2.0f / static_cast<float>(Length));
 
             float real = math::ToFloat(fftResult[k].Real());
             float imag = math::ToFloat(fftResult[k].Imaginary());
@@ -76,16 +76,18 @@ namespace analysis
     }
 
     template<typename QNumberType, std::size_t Length>
-    typename DiscreteConsineTransform<QNumberType, Length>::VectorReal& DiscreteConsineTransform<QNumberType, Length>::Inverse(VectorReal& input)
+    OPTIMIZE_FOR_SPEED
+        typename DiscreteCosineTransform<QNumberType, Length>::VectorReal&
+        DiscreteCosineTransform<QNumberType, Length>::Inverse(VectorReal& input)
     {
         float sqrtN = math::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 real = math::ToFloat(input[k]) * sqrtN / math::Sqrt(2.0f);
+            float imag = -math::ToFloat(input[Length - k]) * sqrtN / math::Sqrt(2.0f);
 
             float angle = static_cast<float>(k) * std::numbers::pi_v<float> / (2.0f * static_cast<float>(Length));
             float cosine = math::Cos(angle);
@@ -106,8 +108,8 @@ namespace analysis
     }
 
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
-    extern template class DiscreteConsineTransform<float, 8>;
-    extern template class DiscreteConsineTransform<math::Q15, 8>;
-    extern template class DiscreteConsineTransform<math::Q31, 8>;
+    extern template class DiscreteCosineTransform<float, 8>;
+    extern template class DiscreteCosineTransform<math::Q15, 8>;
+    extern template class DiscreteCosineTransform<math::Q31, 8>;
 #endif
 }

numerical/analysis/PowerDensitySpectrum.hpp

@@ -85,7 +85,10 @@ namespace analysis
             auto& spectrum = fft.Forward(segment);
 
             for (std::size_t k = 0; k <= SegmentSize / 2; ++k)
-                y[k] += QNumberType(math::ToFloat(MagnitudeSquared(spectrum[k])) / static_cast<float>(SegmentSize));
+            {
+                float factor = (k == 0 || k == SegmentSize / 2) ? 1.0f : 2.0f;
+                y[k] += QNumberType(math::ToFloat(MagnitudeSquared(spectrum[k])) * factor / static_cast<float>(SegmentSize));
+            }
 
             ++segmentCount;
         }

numerical/analysis/test/PowerDensitySpectrumTestSupport.hpp

@@ -70,4 +70,42 @@ namespace analysis::test
         typename VectorComplex::template WithMaxSize<Length> result;
         typename VectorReal::template WithMaxSize<Length> timeResult;
     };
+
+    template<typename QNumberType, std::size_t Length>
+    class FftStubInteriorBin
+        : public analysis::FastFourierTransform<QNumberType>
+    {
+    public:
+        using VectorComplex = typename analysis::FastFourierTransform<QNumberType>::VectorComplex;
+        using VectorReal = typename analysis::FastFourierTransform<QNumberType>::VectorReal;
+
+        explicit FftStubInteriorBin(analysis::TwiddleFactors<QNumberType, Length / 2>&)
+        {}
+
+        VectorComplex& Forward(VectorReal& input) override
+        {
+            result.clear();
+
+            for (std::size_t i = 0; i < Length; ++i)
+            {
+                if (i == 1)
+                    result.push_back(math::Complex<QNumberType>(QNumberType(0.5f), QNumberType(0.0f)));
+                else
+                    result.push_back(math::Complex<QNumberType>(QNumberType(0.0f), QNumberType(0.0f)));
+            }
+            return result;
+        }
+
+        VectorReal& Inverse(VectorComplex& input) override
+        {
+            timeResult.clear();
+            for (std::size_t i = 0; i < Length; ++i)
+                timeResult.push_back(QNumberType(0.0f));
+            return timeResult;
+        }
+
+    private:
+        typename VectorComplex::template WithMaxSize<Length> result;
+        typename VectorReal::template WithMaxSize<Length> timeResult;
+    };
 }

numerical/analysis/test/TestDiscreteCosineTransform.cpp

@@ -49,7 +49,7 @@ namespace
 
         ConcreteTwiddleFactors<float, Length / 2> twiddle;
         analysis::FastFourierTransformRadix2Impl<float, Length> fft{ twiddle };
-        analysis::DiscreteConsineTransform<float, Length> dct{ fft };
+        analysis::DiscreteCosineTransform<float, Length> dct{ fft };
 
         RealBuf signal;
 
@@ -71,7 +71,7 @@ namespace
         using ComplexBuf = VectorComplex::WithMaxSize<Length>;
 
         ::testing::StrictMock<MockFft> mockFft;
-        analysis::DiscreteConsineTransform<float, Length> dct{ mockFft };
+        analysis::DiscreteCosineTransform<float, Length> dct{ mockFft };
 
         RealBuf signal;
         ComplexBuf fftOut;
@@ -112,7 +112,7 @@ TEST_F(TestDiscreteCosineTransform, forward_impulse_at_origin_matches_closed_for
     EXPECT_NEAR(result[0], 1.0f / sqrtN, math::Tolerance<float>());
     for (std::size_t k = 1; k < Length; ++k)
     {
-        float ref{ 2.0f * std::cos(static_cast<float>(k) * std::numbers::pi_v<float> / (2.0f * static_cast<float>(Length))) / sqrtN };
+        float ref{ std::sqrt(2.0f / static_cast<float>(Length)) * std::cos(static_cast<float>(k) * std::numbers::pi_v<float> / (2.0f * static_cast<float>(Length))) };
         EXPECT_NEAR(result[k], ref, math::Tolerance<float>());
     }
 }
@@ -136,7 +136,7 @@ TEST_F(TestDiscreteCosineTransform, forward_matches_direct_dct_ii_definition)
         for (std::size_t n = 0; n < Length; ++n)
             acc += x[n] * std::cos(std::numbers::pi_v<float> * (2.0f * static_cast<float>(n) + 1.0f) * static_cast<float>(k) / (2.0f * static_cast<float>(Length)));
 
-        float ref{ 2.0f / std::sqrt(static_cast<float>(Length)) * acc };
+        float ref{ std::sqrt(2.0f / static_cast<float>(Length)) * acc };
         EXPECT_NEAR(result[k], ref, math::Tolerance<float>());
     }
 }
@@ -224,6 +224,25 @@ TEST_F(TestDiscreteCosineTransform, inverse_pure_dc_spectrum_gives_constant_sign
         EXPECT_NEAR(result[n], ref, math::Tolerance<float>());
 }
 
+TEST_F(TestDiscreteCosineTransform, parseval_identity_holds)
+{
+    constexpr std::array<float, Length> x{ 3.0f, -1.0f, 4.0f, 1.0f, -5.0f, 9.0f, -2.0f, 6.0f };
+    for (std::size_t n = 0; n < Length; ++n)
+        signal[n] = x[n];
+
+    auto& X = dct.Forward(signal);
+
+    float energyTime{ 0.0f };
+    for (std::size_t n = 0; n < Length; ++n)
+        energyTime += x[n] * x[n];
+
+    float energyFreq{ 0.0f };
+    for (std::size_t k = 0; k < Length; ++k)
+        energyFreq += X[k] * X[k];
+
+    EXPECT_NEAR(energyFreq, energyTime, 1.0f);
+}
+
 TEST_F(TestDiscreteCosineTransformMockFft, forward_passes_full_length_buffer_to_fft)
 {
     for (std::size_t i = 0; i < Length; ++i)

numerical/analysis/test/TestPowerDensitySpectrum.cpp

@@ -176,6 +176,34 @@ TYPED_TEST(TestPowerSpectralDensityZeroOverlap, output_size_equals_half_segment_
     EXPECT_EQ(result.size(), TestFixture::length / 2 + 1);
 }
 
+TYPED_TEST(TestPowerSpectralDensity, interior_bin_is_doubled_relative_to_dc_for_equal_magnitude)
+{
+    if constexpr (!std::is_floating_point_v<TypeParam>)
+        GTEST_SKIP();
+
+    using FftInterior = analysis::test::FftStubInteriorBin<TypeParam, TestFixture::length>;
+    using Twiddle = analysis::test::TwiddleFactorsStub<TypeParam, TestFixture::length / 2>;
+    using PsdInterior = analysis::PowerSpectralDensity<TypeParam, TestFixture::length, FftInterior, Twiddle, TestFixture::overlap>;
+
+    analysis::test::WindowStub<TypeParam> win;
+    TypeParam samplingTime = TypeParam(1.0f / 48000.0f);
+
+    PsdInterior psdInterior(win, samplingTime);
+
+    typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize<TestFixture::length> input;
+    for (std::size_t i = 0; i < this->length; ++i)
+        input.push_back(TypeParam(0.5f));
+
+    auto& dcResult = this->powerDensitySpectrum->Calculate(input);
+    float dcBin = math::ToFloat(dcResult[0]);
+
+    auto& interiorResult = psdInterior.Calculate(input);
+    float interiorBin1 = math::ToFloat(interiorResult[1]);
+
+    ASSERT_GT(dcBin, 0.0f);
+    EXPECT_NEAR(interiorBin1 / dcBin, 2.0f, 1e-3f);
+}
+
 TYPED_TEST(TestPowerSpectralDensityZeroOverlap, zero_overlap_step_equals_segment_size)
 {
     if constexpr (!std::is_floating_point_v<TypeParam>)

numerical/control_analysis/TransferFunctionStateSpace.hpp

@@ -150,7 +150,7 @@ namespace control_analysis
         for (std::size_t k{ 1 }; k <= n; ++k)
         {
             auto CB{ sys.C * N * sys.B };
-            numCoeffs[k] = CB.at(0, 0);
+            numCoeffs[k] = CB.at(0, 0) + sys.D.at(0, 0) * charPoly[k];
             N = sys.A * N;
             for (std::size_t i{ 0 }; i < n; ++i)
                 N.at(i, i) += charPoly[k];