From d7ec2f50e719ede5f50949ea113579e91e900256 Mon Sep 17 00:00:00 2001 From: gfs Date: Fri, 14 Aug 2026 20:09:14 +0000 Subject: [PATCH 1/2] fix(analysis,control_analysis,math): resolve issues 251, 252, 253 and add coverage MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - #251: rename DiscreteConsineTransform to DiscreteCosineTransform; fix forward scale from 2/√N to √(2/N) and inverse scale from √N/2 to √N/√2 to produce the orthonormal DCT-II; update doc and two reference-value tests; add Parseval identity test; add OPTIMIZE_FOR_SPEED to Inverse - #252: multiply interior bins (1..N/2-1) by 2 in PowerSpectralDensity::Calculate to correct 3 dB under-read in one-sided PSD; add FftStubInteriorBin and ratio test verifying the doubling factor - #253: add missing D·charPoly[k] term in TransferFunctionStateSpace::ToTransferFunction so systems with nonzero feedthrough return correct numerator coefficients; add round-trip and DC-gain tests for improper transfer functions - Add TestTolerance.cpp with five TEST_F cases covering float, Q31, Q15 - Add seven TYPED_TEST cases to TestQNumber covering zero construction, sign rules, double negation, and compound-divide-by-zero death - Add six TEST_F cases to TestFrequencyResponse covering size, phase, and DC-gain properties across Biquad, Highpass, FirstOrderIir, PureDelay, and Lowpass fixtures Co-Authored-By: Claude Sonnet 4.6 --- doc/analysis/DiscreteCosineTransform.md | 36 +++++----- .../analysis/DiscreteCosineTransform.cpp | 6 +- .../analysis/DiscreteCosineTransform.hpp | 30 ++++---- numerical/analysis/PowerDensitySpectrum.hpp | 5 +- .../test/PowerDensitySpectrumTestSupport.hpp | 38 ++++++++++ .../test/TestDiscreteCosineTransform.cpp | 27 ++++++-- .../test/TestPowerDensitySpectrum.cpp | 28 ++++++++ .../TransferFunctionStateSpace.hpp | 2 +- .../test/TestFrequencyResponse.cpp | 69 ++++++++++++++++++- .../test/TestTransferFunctionStateSpace.cpp | 29 ++++++++ numerical/math/test/CMakeLists.txt | 1 + numerical/math/test/TestQNumber.cpp | 66 ++++++++++++++++++ numerical/math/test/TestTolerance.cpp | 34 +++++++++ 13 files changed, 330 insertions(+), 41 deletions(-) create mode 100644 numerical/math/test/TestTolerance.cpp diff --git a/doc/analysis/DiscreteCosineTransform.md b/doc/analysis/DiscreteCosineTransform.md index 99c6b2a6..c3ba9071 100644 --- a/doc/analysis/DiscreteCosineTransform.md +++ b/doc/analysis/DiscreteCosineTransform.md @@ -10,13 +10,15 @@ 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 @@ -24,7 +26,9 @@ 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. diff --git a/numerical/analysis/DiscreteCosineTransform.cpp b/numerical/analysis/DiscreteCosineTransform.cpp index a0e2998f..fd9466fc 100644 --- a/numerical/analysis/DiscreteCosineTransform.cpp +++ b/numerical/analysis/DiscreteCosineTransform.cpp @@ -3,7 +3,7 @@ namespace analysis { - template class DiscreteConsineTransform; - template class DiscreteConsineTransform; - template class DiscreteConsineTransform; + template class DiscreteCosineTransform; + template class DiscreteCosineTransform; + template class DiscreteCosineTransform; } diff --git a/numerical/analysis/DiscreteCosineTransform.hpp b/numerical/analysis/DiscreteCosineTransform.hpp index 15f53a97..3d13d875 100644 --- a/numerical/analysis/DiscreteCosineTransform.hpp +++ b/numerical/analysis/DiscreteCosineTransform.hpp @@ -14,17 +14,17 @@ namespace analysis { template - 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::value || std::is_floating_point_v, - "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::VectorReal; using VectorComplex = typename FastFourierTransform::VectorComplex; - explicit DiscreteConsineTransform(FastFourierTransform& fft); + explicit DiscreteCosineTransform(FastFourierTransform& fft); VectorReal& Forward(VectorReal& input); VectorReal& Inverse(VectorReal& input); @@ -38,7 +38,7 @@ namespace analysis // Implementation // template - DiscreteConsineTransform::DiscreteConsineTransform(FastFourierTransform& fft) + DiscreteCosineTransform::DiscreteCosineTransform(FastFourierTransform& fft) : fft(fft) { output.resize(Length); @@ -48,8 +48,8 @@ namespace analysis template OPTIMIZE_FOR_SPEED - typename DiscreteConsineTransform::VectorReal& - DiscreteConsineTransform::Forward(VectorReal& input) + typename DiscreteCosineTransform::VectorReal& + DiscreteCosineTransform::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(k) * std::numbers::pi_v / (2.0f * Length); - float scale = 2.0f / math::Sqrt(static_cast(Length)); + float scale = math::Sqrt(2.0f / static_cast(Length)); float real = math::ToFloat(fftResult[k].Real()); float imag = math::ToFloat(fftResult[k].Imaginary()); @@ -76,7 +76,9 @@ namespace analysis } template - typename DiscreteConsineTransform::VectorReal& DiscreteConsineTransform::Inverse(VectorReal& input) + OPTIMIZE_FOR_SPEED + typename DiscreteCosineTransform::VectorReal& + DiscreteCosineTransform::Inverse(VectorReal& input) { float sqrtN = math::Sqrt(static_cast(Length)); @@ -84,8 +86,8 @@ namespace analysis 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(k) * std::numbers::pi_v / (2.0f * static_cast(Length)); float cosine = math::Cos(angle); @@ -106,8 +108,8 @@ namespace analysis } #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD - extern template class DiscreteConsineTransform; - extern template class DiscreteConsineTransform; - extern template class DiscreteConsineTransform; + extern template class DiscreteCosineTransform; + extern template class DiscreteCosineTransform; + extern template class DiscreteCosineTransform; #endif } diff --git a/numerical/analysis/PowerDensitySpectrum.hpp b/numerical/analysis/PowerDensitySpectrum.hpp index f7e66908..8901c5bd 100644 --- a/numerical/analysis/PowerDensitySpectrum.hpp +++ b/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(SegmentSize)); + { + float factor = (k == 0 || k == SegmentSize / 2) ? 1.0f : 2.0f; + y[k] += QNumberType(math::ToFloat(MagnitudeSquared(spectrum[k])) * factor / static_cast(SegmentSize)); + } ++segmentCount; } diff --git a/numerical/analysis/test/PowerDensitySpectrumTestSupport.hpp b/numerical/analysis/test/PowerDensitySpectrumTestSupport.hpp index 8f9a9b41..8e8f3e26 100644 --- a/numerical/analysis/test/PowerDensitySpectrumTestSupport.hpp +++ b/numerical/analysis/test/PowerDensitySpectrumTestSupport.hpp @@ -70,4 +70,42 @@ namespace analysis::test typename VectorComplex::template WithMaxSize result; typename VectorReal::template WithMaxSize timeResult; }; + + template + class FftStubInteriorBin + : public analysis::FastFourierTransform + { + public: + using VectorComplex = typename analysis::FastFourierTransform::VectorComplex; + using VectorReal = typename analysis::FastFourierTransform::VectorReal; + + explicit FftStubInteriorBin(analysis::TwiddleFactors&) + {} + + 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(0.5f), QNumberType(0.0f))); + else + result.push_back(math::Complex(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 result; + typename VectorReal::template WithMaxSize timeResult; + }; } diff --git a/numerical/analysis/test/TestDiscreteCosineTransform.cpp b/numerical/analysis/test/TestDiscreteCosineTransform.cpp index 15cf794d..ed999292 100644 --- a/numerical/analysis/test/TestDiscreteCosineTransform.cpp +++ b/numerical/analysis/test/TestDiscreteCosineTransform.cpp @@ -49,7 +49,7 @@ namespace ConcreteTwiddleFactors twiddle; analysis::FastFourierTransformRadix2Impl fft{ twiddle }; - analysis::DiscreteConsineTransform dct{ fft }; + analysis::DiscreteCosineTransform dct{ fft }; RealBuf signal; @@ -71,7 +71,7 @@ namespace using ComplexBuf = VectorComplex::WithMaxSize; ::testing::StrictMock mockFft; - analysis::DiscreteConsineTransform dct{ mockFft }; + analysis::DiscreteCosineTransform 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()); for (std::size_t k = 1; k < Length; ++k) { - float ref{ 2.0f * std::cos(static_cast(k) * std::numbers::pi_v / (2.0f * static_cast(Length))) / sqrtN }; + float ref{ std::sqrt(2.0f / static_cast(Length)) * std::cos(static_cast(k) * std::numbers::pi_v / (2.0f * static_cast(Length))) }; EXPECT_NEAR(result[k], ref, math::Tolerance()); } } @@ -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 * (2.0f * static_cast(n) + 1.0f) * static_cast(k) / (2.0f * static_cast(Length))); - float ref{ 2.0f / std::sqrt(static_cast(Length)) * acc }; + float ref{ std::sqrt(2.0f / static_cast(Length)) * acc }; EXPECT_NEAR(result[k], ref, math::Tolerance()); } } @@ -224,6 +224,25 @@ TEST_F(TestDiscreteCosineTransform, inverse_pure_dc_spectrum_gives_constant_sign EXPECT_NEAR(result[n], ref, math::Tolerance()); } +TEST_F(TestDiscreteCosineTransform, parseval_identity_holds) +{ + constexpr std::array 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) diff --git a/numerical/analysis/test/TestPowerDensitySpectrum.cpp b/numerical/analysis/test/TestPowerDensitySpectrum.cpp index 906797e9..0cc810be 100644 --- a/numerical/analysis/test/TestPowerDensitySpectrum.cpp +++ b/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) + GTEST_SKIP(); + + using FftInterior = analysis::test::FftStubInteriorBin; + using Twiddle = analysis::test::TwiddleFactorsStub; + using PsdInterior = analysis::PowerSpectralDensity; + + analysis::test::WindowStub win; + TypeParam samplingTime = TypeParam(1.0f / 48000.0f); + + PsdInterior psdInterior(win, samplingTime); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize 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) diff --git a/numerical/control_analysis/TransferFunctionStateSpace.hpp b/numerical/control_analysis/TransferFunctionStateSpace.hpp index afd84afd..ee653838 100644 --- a/numerical/control_analysis/TransferFunctionStateSpace.hpp +++ b/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]; diff --git a/numerical/control_analysis/test/TestFrequencyResponse.cpp b/numerical/control_analysis/test/TestFrequencyResponse.cpp index 4f156f71..4778903e 100644 --- a/numerical/control_analysis/test/TestFrequencyResponse.cpp +++ b/numerical/control_analysis/test/TestFrequencyResponse.cpp @@ -3,13 +3,11 @@ #include #include #include -#include #include namespace { static constexpr float kSampleFrequency = 1000.0f; - static constexpr float kPi = std::numbers::pi_v; class TestFrequencyResponseUnity : public ::testing::Test { @@ -332,3 +330,70 @@ TEST_F(TestFrequencyResponseUnity, zero_denominator_coefficients_produce_finite_ EXPECT_FALSE(std::isinf(m)); } } + +TEST_F(TestFrequencyResponseUnity, magnitude_output_size_matches_points_for_128_points) +{ + std::array b{ 1.0f }; + std::array a{ 1.0f }; + control_analysis::FrequencyResponse fr128{ b, a, kSampleFrequency }; + + auto [frequencies, magnitudes, phases] = fr128.Calculate(); + + EXPECT_EQ(magnitudes.size(), 128u); + EXPECT_EQ(phases.size(), 128u); + EXPECT_EQ(frequencies.size(), 128u); +} + +TEST_F(TestFrequencyResponseHighpass, phase_is_positive_at_low_frequencies) +{ + auto [frequencies, magnitudes, phases] = freqResponse.Calculate(); + + ASSERT_FALSE(phases.empty()); + EXPECT_GT(phases.front(), 0.0f); +} + +TEST_F(TestFrequencyResponseBiquad, dc_gain_is_finite_and_nonnan) +{ + auto [frequencies, magnitudes, phases] = freqResponse.Calculate(); + + ASSERT_FALSE(magnitudes.empty()); + EXPECT_FALSE(std::isnan(magnitudes.front())); + EXPECT_FALSE(std::isinf(magnitudes.front())); +} + +TEST_F(TestFrequencyResponsePureDelay, phase_at_nyquist_is_minus_180_degrees) +{ + auto [frequencies, magnitudes, phases] = freqResponse.Calculate(); + + ASSERT_FALSE(phases.empty()); + EXPECT_NEAR(std::abs(phases.back()), 180.0f, 5.0f); +} + +TEST_F(TestFrequencyResponseFirstOrderIir, phase_is_negative_at_quarter_nyquist) +{ + auto [frequencies, magnitudes, phases] = freqResponse.Calculate(); + + float targetFreq = kSampleFrequency / 4.0f; + float bestPhase = 0.0f; + float minDiff = std::numeric_limits::max(); + + for (std::size_t i = 0; i < frequencies.size(); ++i) + { + float diff = std::abs(frequencies[i] - targetFreq); + if (diff < minDiff) + { + minDiff = diff; + bestPhase = phases[i]; + } + } + + EXPECT_NEAR(bestPhase, -26.57f, 2.0f); +} + +TEST_F(TestFrequencyResponseLowpass, lowpass_dc_gain_is_zero_db) +{ + auto [frequencies, magnitudes, phases] = freqResponse.Calculate(); + + ASSERT_FALSE(magnitudes.empty()); + EXPECT_NEAR(magnitudes.front(), 0.0f, 0.1f); +} diff --git a/numerical/control_analysis/test/TestTransferFunctionStateSpace.cpp b/numerical/control_analysis/test/TestTransferFunctionStateSpace.cpp index 4e0e4261..55b6c63a 100644 --- a/numerical/control_analysis/test/TestTransferFunctionStateSpace.cpp +++ b/numerical/control_analysis/test/TestTransferFunctionStateSpace.cpp @@ -166,3 +166,32 @@ TEST_F(TestTransferFunctionStateSpace, non_monic_round_trip_recovers_monic_tf) EXPECT_NEAR(tf2.denominator[1], 3.0f, math::Tolerance()); EXPECT_NEAR(tf2.denominator[2], 2.0f, math::Tolerance()); } + +TEST_F(TestTransferFunctionStateSpace, round_trip_with_nonzero_feedthrough_recovers_numerator) +{ + TF2 improper{}; + improper.denominator = { 1.0f, 3.0f, 2.0f }; + improper.numerator = { 2.0f, 7.0f, 9.0f }; + + auto tf2{ TFSS::ToTransferFunction(TFSS::ToControllableCanonical(improper)) }; + + EXPECT_NEAR(tf2.numerator[0], 2.0f, 1e-4f); + EXPECT_NEAR(tf2.numerator[1], 7.0f, 1e-4f); + EXPECT_NEAR(tf2.numerator[2], 9.0f, 1e-4f); + EXPECT_NEAR(tf2.denominator[0], 1.0f, 1e-4f); + EXPECT_NEAR(tf2.denominator[1], 3.0f, 1e-4f); + EXPECT_NEAR(tf2.denominator[2], 2.0f, 1e-4f); +} + +TEST_F(TestTransferFunctionStateSpace, dc_gain_correct_for_system_with_feedthrough) +{ + TF2 improper{}; + improper.denominator = { 1.0f, 1.0f, 0.0f }; + improper.numerator = { 2.0f, 3.0f, 0.0f }; + + auto sys{ TFSS::ToControllableCanonical(improper) }; + auto tf2{ TFSS::ToTransferFunction(sys) }; + + EXPECT_NEAR(tf2.numerator[0], 2.0f, 1e-4f); + EXPECT_NEAR(tf2.numerator[1], 3.0f, 1e-4f); +} diff --git a/numerical/math/test/CMakeLists.txt b/numerical/math/test/CMakeLists.txt index bb07325a..c68a8310 100644 --- a/numerical/math/test/CMakeLists.txt +++ b/numerical/math/test/CMakeLists.txt @@ -23,6 +23,7 @@ target_sources(numerical.math_test PRIVATE TestMatrixNorms.cpp TestMatrixOperations.cpp TestQNumber.cpp + TestTolerance.cpp TestMatrix.cpp TestQuaternion.cpp TestRecursiveBuffer.cpp diff --git a/numerical/math/test/TestQNumber.cpp b/numerical/math/test/TestQNumber.cpp index ed921ba0..a0b59e42 100644 --- a/numerical/math/test/TestQNumber.cpp +++ b/numerical/math/test/TestQNumber.cpp @@ -353,3 +353,69 @@ TEST_F(QNumberUtilTest, MinMaxLowest_Q15) EXPECT_NEAR(math::Max(), 0.9999f, 1e-3f); EXPECT_NEAR(math::Lowest(), -0.9999f, 1e-3f); } + +TYPED_TEST(QNumberTest, FloatConstructorZeroIsZero) +{ + TypeParam num(0.0f); + + EXPECT_EQ(num.RawValue(), 0); + EXPECT_FLOAT_EQ(num.ToFloat(), 0.0f); +} + +TYPED_TEST(QNumberTest, NegativeTimesNegativeIsPositive) +{ + TypeParam a(-0.25f); + TypeParam b(-0.25f); + + TypeParam result = a * b; + + EXPECT_GT(result.ToFloat(), 0.0f); + EXPECT_NEAR(result.ToFloat(), 0.0625f, math::Tolerance()); +} + +TYPED_TEST(QNumberTest, NegativeDividedByNegativeIsPositive) +{ + TypeParam a(-0.20f); + TypeParam b(-0.40f); + + TypeParam result = a / b; + + EXPECT_NEAR(result.ToFloat(), 0.50f, math::Tolerance()); +} + +TYPED_TEST(QNumberTest, UnaryNegationOfZeroIsZero) +{ + TypeParam zero; + + TypeParam result = -zero; + + EXPECT_EQ(result.RawValue(), 0); +} + +TYPED_TEST(QNumberTest, SubtractLargerFromSmallerIsNegative) +{ + TypeParam a(0.10f); + TypeParam b(0.30f); + + TypeParam result = a - b; + + EXPECT_LT(result.ToFloat(), 0.0f); + EXPECT_NEAR(result.ToFloat(), -0.20f, math::Tolerance()); +} + +TYPED_TEST(QNumberTest, DoubleNegationIsIdentity) +{ + TypeParam a(0.20f); + + TypeParam result = -(-a); + + EXPECT_NEAR(result.ToFloat(), a.ToFloat(), math::Tolerance()); +} + +TYPED_TEST(QNumberTest, CompoundDivideByZeroDies) +{ + TypeParam a(0.10f); + TypeParam zero(0.0f); + + EXPECT_DEATH_IF_SUPPORTED({ a /= zero; }, ""); // NOLINT +} diff --git a/numerical/math/test/TestTolerance.cpp b/numerical/math/test/TestTolerance.cpp new file mode 100644 index 00000000..b463ae9c --- /dev/null +++ b/numerical/math/test/TestTolerance.cpp @@ -0,0 +1,34 @@ +#include "numerical/math/QNumber.hpp" +#include "numerical/math/Tolerance.hpp" +#include + +namespace +{ + class TestTolerance : public ::testing::Test + {}; +} + +TEST_F(TestTolerance, tolerance_for_float_returns_1e_minus_3) +{ + EXPECT_FLOAT_EQ(math::Tolerance(), 1e-3f); +} + +TEST_F(TestTolerance, tolerance_for_q31_returns_1e_minus_3) +{ + EXPECT_FLOAT_EQ(math::Tolerance(), 1e-3f); +} + +TEST_F(TestTolerance, tolerance_for_q15_returns_1e_minus_3) +{ + EXPECT_FLOAT_EQ(math::Tolerance(), 1e-3f); +} + +TEST_F(TestTolerance, tolerance_is_positive) +{ + EXPECT_GT(math::Tolerance(), 0.0f); +} + +TEST_F(TestTolerance, tolerance_is_less_than_one) +{ + EXPECT_LT(math::Tolerance(), 1.0f); +} From c7d28c670ecbc70bfab33b030fb2d27c91dbed6d Mon Sep 17 00:00:00 2001 From: gfs Date: Fri, 14 Aug 2026 22:14:50 +0200 Subject: [PATCH 2/2] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- doc/analysis/DiscreteCosineTransform.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/analysis/DiscreteCosineTransform.md b/doc/analysis/DiscreteCosineTransform.md index c3ba9071..3b50264c 100644 --- a/doc/analysis/DiscreteCosineTransform.md +++ b/doc/analysis/DiscreteCosineTransform.md @@ -58,11 +58,11 @@ $$Y = \text{FFT}([1, 3, 4, 2]) = [10,\; -3-j,\; 0,\; -3+j]$$ **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]$ | -|-----|----------------|-------------------------|---------------------------------| -| 0 | — | — | $10 / \sqrt{4} = 5$ | +| $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$ | +| 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 [5,\; -2.23,\; 0,\; -0.16]$