From e3f738189c5134cfaf2a6f04bfd924cef183702d Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 11:23:14 +0000 Subject: [PATCH 1/9] improve DCT --- .../analysis/DiscreteCosineTransform.hpp | 33 +- .../test/TestDiscreteCosineTransform.cpp | 286 ++++++++++++++---- 2 files changed, 248 insertions(+), 71 deletions(-) diff --git a/numerical/analysis/DiscreteCosineTransform.hpp b/numerical/analysis/DiscreteCosineTransform.hpp index 72a5d022..8dfdae73 100644 --- a/numerical/analysis/DiscreteCosineTransform.hpp +++ b/numerical/analysis/DiscreteCosineTransform.hpp @@ -31,6 +31,7 @@ namespace analysis private: FastFourierTransform& fft; typename infra::BoundedVector::template WithMaxSize output; + typename infra::BoundedVector::template WithMaxSize reordered; typename VectorComplex::template WithMaxSize 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::VectorReal& DiscreteConsineTransform::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(Length))); @@ -70,18 +78,31 @@ namespace analysis template typename DiscreteConsineTransform::VectorReal& DiscreteConsineTransform::Inverse(VectorReal& input) { - complexBuffer[0] = math::Complex{ QNumberType(math::ToFloat(input[0]) * std::sqrt(static_cast(Length))), QNumberType(0.0f) }; + float sqrtN = std::sqrt(static_cast(Length)); + + complexBuffer[0] = math::Complex{ 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(k) * std::numbers::pi_v / (2.0f * static_cast(Length)); - float scale = std::sqrt(static_cast(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(value * std::cos(angle)), QNumberType(value * std::sin(angle)) }; + complexBuffer[k] = math::Complex{ 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 diff --git a/numerical/analysis/test/TestDiscreteCosineTransform.cpp b/numerical/analysis/test/TestDiscreteCosineTransform.cpp index 2adc34b6..a93df39b 100644 --- a/numerical/analysis/test/TestDiscreteCosineTransform.cpp +++ b/numerical/analysis/test/TestDiscreteCosineTransform.cpp @@ -1,116 +1,272 @@ #include "numerical/analysis/DiscreteCosineTransform.hpp" +#include "numerical/analysis/FastFourierTransformRadix2Impl.hpp" #include "numerical/math/QNumber.hpp" +#include "numerical/math/Tolerance.hpp" #include "gmock/gmock.h" +#include +#include +#include namespace { - template - class MockFft - : public analysis::FastFourierTransform + template + class ConcreteTwiddleFactors : public analysis::TwiddleFactors { public: - MOCK_METHOD(typename analysis::FastFourierTransform::VectorComplex&, Forward, (typename analysis::FastFourierTransform::VectorReal & input), (override)); - MOCK_METHOD(typename analysis::FastFourierTransform::VectorReal&, Inverse, (typename analysis::FastFourierTransform::VectorComplex & input), (override)); + ConcreteTwiddleFactors() + { + for (std::size_t k = 0; k < HalfLen; ++k) + { + float angle{ -2.0f * std::numbers::pi_v * static_cast(k) / static_cast(2 * HalfLen) }; + factors[k] = math::Complex{ T(std::cos(angle)), T(std::sin(angle)) }; + } + } + + math::Complex& operator[](std::size_t n) override + { + return factors[n]; + } + + private: + std::array, HalfLen> factors; + }; + + class MockFft : public analysis::FastFourierTransform + { + public: + MOCK_METHOD(VectorComplex&, Forward, (VectorReal& input), (override)); + MOCK_METHOD(VectorReal&, Inverse, (VectorComplex& input), (override)); }; - template - class TestDiscreteConsineTransform - : public ::testing::Test + class TestDiscreteCosineTransform : public ::testing::Test { public: static constexpr std::size_t Length = 8; - using VectorReal = typename analysis::FastFourierTransform::VectorReal; - using VectorComplex = typename analysis::FastFourierTransform::VectorComplex; - ::testing::StrictMock> mockFft; - std::optional> dct; - typename VectorReal::template WithMaxSize input; - typename VectorReal::template WithMaxSize output; - typename VectorComplex::template WithMaxSize fftOutput; + using VectorReal = analysis::FastFourierTransform::VectorReal; + using VectorComplex = analysis::FastFourierTransform::VectorComplex; + using RealBuf = VectorReal::WithMaxSize; + + ConcreteTwiddleFactors twiddle; + analysis::FastFourierTransformRadix2Impl fft{ twiddle }; + analysis::DiscreteConsineTransform dct{ fft }; + + RealBuf signal; void SetUp() override { - dct.emplace(mockFft); - input.clear(); - input.resize(Length); + signal.clear(); + signal.resize(Length); } + }; - typename VectorComplex::template WithMaxSize& EmptyForwardResult() - { - fftOutput.clear(); - fftOutput.resize(fftOutput.max_size()); - return fftOutput; - } + class TestDiscreteCosineTransformMockFft : public ::testing::Test + { + public: + static constexpr std::size_t Length = 8; + + using VectorReal = analysis::FastFourierTransform::VectorReal; + using VectorComplex = analysis::FastFourierTransform::VectorComplex; + using RealBuf = VectorReal::WithMaxSize; + using ComplexBuf = VectorComplex::WithMaxSize; + + ::testing::StrictMock mockFft; + analysis::DiscreteConsineTransform dct{ mockFft }; - typename VectorReal::template WithMaxSize& EmptyInverseResult() + RealBuf signal; + ComplexBuf fftOut; + RealBuf ifftOut; + + void SetUp() override { - output.clear(); - output.resize(output.max_size()); - return output; + signal.clear(); + signal.resize(Length); + + fftOut.clear(); + fftOut.resize(Length); + + ifftOut.clear(); + ifftOut.resize(Length); } }; +} + +TEST_F(TestDiscreteCosineTransform, forward_dc_signal_concentrates_energy_in_bin_zero) +{ + std::fill(signal.begin(), signal.end(), 1.0f); + + auto& result = dct.Forward(signal); - using TestedTypes = ::testing::Types; - TYPED_TEST_SUITE(TestDiscreteConsineTransform, TestedTypes); + EXPECT_NEAR(result[0], std::sqrt(static_cast(Length)), math::Tolerance()); + for (std::size_t k = 1; k < Length; ++k) + EXPECT_NEAR(result[k], 0.0f, math::Tolerance()); } -TYPED_TEST(TestDiscreteConsineTransform, forward_transform_calls_fft_with_preprocessed_data) +TEST_F(TestDiscreteCosineTransform, forward_impulse_at_origin_matches_closed_form) { - for (size_t i = 0; i < this->input.size(); ++i) - this->input[i] = TypeParam(0.1f); + signal[0] = 1.0f; - EXPECT_CALL(this->mockFft, Forward(::testing::_)) - .WillOnce(::testing::Invoke([this](typename analysis::FastFourierTransform::VectorReal& input) -> typename analysis::FastFourierTransform::VectorComplex& - { - for (size_t i = 0; i < input.size() / 2; ++i) - EXPECT_EQ(math::ToFloat(input[i]), math::ToFloat(this->input[2 * i])); + auto& result = dct.Forward(signal); - return this->EmptyForwardResult(); - })); + constexpr float sqrtN{ 2.82842712f }; + 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 }; + EXPECT_NEAR(result[k], ref, math::Tolerance()); + } +} + +TEST_F(TestDiscreteCosineTransform, forward_matches_direct_dct_ii_definition) +{ + 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& result = dct.Forward(signal); - this->dct->Forward(this->input); + float sum{ 0.0f }; + for (std::size_t n = 0; n < Length; ++n) + sum += x[n]; + EXPECT_NEAR(result[0], sum / std::sqrt(static_cast(Length)), math::Tolerance()); + + for (std::size_t k = 1; k < Length; ++k) + { + float acc{ 0.0f }; + 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 }; + EXPECT_NEAR(result[k], ref, math::Tolerance()); + } } -TYPED_TEST(TestDiscreteConsineTransform, inverse_transform_calls_fft_with_preprocessed_data) +TEST_F(TestDiscreteCosineTransform, inverse_recovers_forward_input) { - for (size_t i = 0; i < this->input.size(); ++i) - this->input[i] = TypeParam(0.1f); + constexpr std::array x{ 0.3f, -0.1f, 0.4f, 0.1f, -0.5f, 0.2f, -0.2f, 0.15f }; + for (std::size_t n = 0; n < Length; ++n) + signal[n] = x[n]; - EXPECT_CALL(this->mockFft, Inverse(::testing::_)) - .WillOnce(::testing::Invoke([this](typename analysis::FastFourierTransform::VectorComplex& input) -> typename analysis::FastFourierTransform::VectorReal& - { - EXPECT_LE(std::abs(math::ToFloat(input[0].Real())), 1.0f); + auto& coeffs = dct.Forward(signal); - return this->EmptyInverseResult(); - })); + RealBuf spectrum; + spectrum.resize(Length); + for (std::size_t k = 0; k < Length; ++k) + spectrum[k] = coeffs[k]; + + auto& recovered = dct.Inverse(spectrum); + + for (std::size_t n = 0; n < Length; ++n) + EXPECT_NEAR(recovered[n], x[n], math::Tolerance()); +} + +TEST_F(TestDiscreteCosineTransform, forward_zero_input_gives_zero_output) +{ + auto& result = dct.Forward(signal); + + for (std::size_t k = 0; k < Length; ++k) + EXPECT_NEAR(result[k], 0.0f, math::Tolerance()); +} + +TEST_F(TestDiscreteCosineTransform, forward_is_linear) +{ + constexpr std::array x1{ 1.0f, 0.0f, -1.0f, 0.0f, 1.0f, 0.0f, -1.0f, 0.0f }; + constexpr std::array x2{ 0.0f, 1.0f, 0.0f, -1.0f, 0.0f, 1.0f, 0.0f, -1.0f }; + constexpr float a{ 0.5f }; + constexpr float b{ -0.25f }; + + VectorReal::WithMaxSize buf1; + VectorReal::WithMaxSize buf2; + VectorReal::WithMaxSize bufComb; + buf1.resize(Length); + buf2.resize(Length); + bufComb.resize(Length); + + for (std::size_t n = 0; n < Length; ++n) + { + buf1[n] = x1[n]; + buf2[n] = x2[n]; + bufComb[n] = a * x1[n] + b * x2[n]; + } + + auto& r1 = dct.Forward(buf1); + std::array dct1{}; + for (std::size_t k = 0; k < Length; ++k) + dct1[k] = r1[k]; - this->dct->Inverse(this->input); + auto& r2 = dct.Forward(buf2); + std::array dct2{}; + for (std::size_t k = 0; k < Length; ++k) + dct2[k] = r2[k]; + + auto& rComb = dct.Forward(bufComb); + + for (std::size_t k = 0; k < Length; ++k) + EXPECT_NEAR(rComb[k], a * dct1[k] + b * dct2[k], math::Tolerance()); +} + +TEST_F(TestDiscreteCosineTransform, inverse_zero_input_gives_zero_output) +{ + auto& result = dct.Inverse(signal); + + for (std::size_t n = 0; n < Length; ++n) + EXPECT_NEAR(result[n], 0.0f, math::Tolerance()); +} + +TEST_F(TestDiscreteCosineTransform, inverse_pure_dc_spectrum_gives_constant_signal) +{ + signal[0] = 1.0f; + + auto& result = dct.Inverse(signal); + + float ref{ 1.0f / std::sqrt(static_cast(Length)) }; + for (std::size_t n = 0; n < Length; ++n) + EXPECT_NEAR(result[n], ref, math::Tolerance()); } -TYPED_TEST(TestDiscreteConsineTransform, input_size_matches_fft_size) +TEST_F(TestDiscreteCosineTransformMockFft, forward_passes_full_length_buffer_to_fft) { - EXPECT_CALL(this->mockFft, Forward(::testing::_)) - .WillOnce(::testing::Invoke([this](typename analysis::FastFourierTransform::VectorReal& input) -> typename analysis::FastFourierTransform::VectorComplex& + for (std::size_t i = 0; i < Length; ++i) + signal[i] = 0.1f; + + EXPECT_CALL(mockFft, Forward(::testing::_)) + .WillOnce(::testing::Invoke([this](VectorReal& in) -> VectorComplex& { - EXPECT_EQ(input.size(), TestDiscreteConsineTransform::Length); - return this->EmptyForwardResult(); + EXPECT_EQ(in.size(), Length); + return fftOut; })); - this->dct->Forward(this->input); + dct.Forward(signal); } -TYPED_TEST(TestDiscreteConsineTransform, scaling_stays_within_fixed_point_range) +TEST_F(TestDiscreteCosineTransformMockFft, inverse_builds_complex_buffer_within_valid_range) { - std::fill(this->input.begin(), this->input.end(), TypeParam(0.9999f)); + for (std::size_t i = 0; i < Length; ++i) + signal[i] = 0.1f; - EXPECT_CALL(this->mockFft, Forward(::testing::_)) - .WillOnce(::testing::Invoke([this](typename analysis::FastFourierTransform::VectorReal& input) -> typename analysis::FastFourierTransform::VectorComplex& + EXPECT_CALL(mockFft, Inverse(::testing::_)) + .WillOnce(::testing::Invoke([this](VectorComplex& in) -> VectorReal& { - for (const auto& value : input) - EXPECT_LE(std::abs(math::ToFloat(value)), 1.0f); + EXPECT_EQ(in.size(), Length); + for (std::size_t k = 0; k < in.size(); ++k) + EXPECT_LE(std::abs(math::ToFloat(in[k].Real())), 1.0f); + return ifftOut; + })); - return this->EmptyForwardResult(); + dct.Inverse(signal); +} + +TEST_F(TestDiscreteCosineTransformMockFft, forward_scaling_stays_within_range_at_full_scale) +{ + std::fill(signal.begin(), signal.end(), 0.9999f); + + EXPECT_CALL(mockFft, Forward(::testing::_)) + .WillOnce(::testing::Invoke([this](VectorReal& in) -> VectorComplex& + { + for (const auto& v : in) + EXPECT_LE(std::abs(math::ToFloat(v)), 1.0f); + return fftOut; })); - this->dct->Forward(this->input); + dct.Forward(signal); } From 3da123cdfef0c38fa21ccedd5ab5a844d41f3c3f Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 11:26:05 +0000 Subject: [PATCH 2/9] improve coverage for fft radix2 --- .../TestFastFourierTransformRadix2Impl.cpp | 177 ++++++++++++++++++ 1 file changed, 177 insertions(+) diff --git a/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp b/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp index 0c57cfa5..046a3fa3 100644 --- a/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp +++ b/numerical/analysis/test/TestFastFourierTransformRadix2Impl.cpp @@ -4,6 +4,7 @@ #include "gtest/gtest.h" #include #include +#include namespace { @@ -59,6 +60,45 @@ namespace using TestedTypes = ::testing::Types; TYPED_TEST_SUITE(TestFastFourierTransform, TestedTypes); + + class RealTwiddleFactors8 : public analysis::TwiddleFactors + { + public: + RealTwiddleFactors8() + { + for (std::size_t k{ 0 }; k < 4; ++k) + { + float angle{ -2.0f * std::numbers::pi_v * static_cast(k) / 8.0f }; + factors[k] = math::Complex{ std::cos(angle), std::sin(angle) }; + } + } + + math::Complex& operator[](std::size_t n) override + { + return factors[n]; + } + + private: + std::array, 4> factors{}; + }; + + class TestFastFourierTransformFloat : public ::testing::Test + { + public: + static constexpr std::size_t Length = 8; + using VectorComplex = analysis::FastFourierTransform::VectorComplex; + using VectorReal = analysis::FastFourierTransform::VectorReal; + + void SetUp() override + { + fft.emplace(twiddleFactors); + } + + RealTwiddleFactors8 twiddleFactors{}; + std::optional> fft; + typename VectorReal::template WithMaxSize timeDomain; + typename VectorComplex::template WithMaxSize frequencyDomain; + }; } TYPED_TEST(TestFastFourierTransform, zero_input_produces_zero_output) @@ -118,3 +158,140 @@ TYPED_TEST(TestFastFourierTransform, nyquist_frequency_detection) static_cast(TestFastFourierTransform::Length) * 0.1f, math::Tolerance()); } + +TEST_F(TestFastFourierTransformFloat, sinusoid_bin1_magnitude_and_phase_match_dft) +{ + constexpr std::size_t k{ 1 }; + timeDomain.resize(Length); + for (std::size_t n{ 0 }; n < Length; ++n) + timeDomain[n] = std::cos(2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(Length)); + + auto& result = fft->Forward(timeDomain); + + float refReal{ 0.0f }; + float refImag{ 0.0f }; + for (std::size_t n{ 0 }; n < Length; ++n) + { + float angle{ -2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(Length) }; + refReal += timeDomain[n] * std::cos(angle); + refImag += timeDomain[n] * std::sin(angle); + } + + EXPECT_NEAR(result[k].Real(), refReal, 1e-3f); + EXPECT_NEAR(result[k].Imaginary(), refImag, 1e-3f); +} + +TEST_F(TestFastFourierTransformFloat, impulse_produces_flat_magnitude_spectrum) +{ + timeDomain.resize(Length, 0.0f); + timeDomain[0] = 1.0f; + + auto& result = fft->Forward(timeDomain); + + for (std::size_t k{ 0 }; k < Length; ++k) + { + float mag{ std::sqrt(result[k].Real() * result[k].Real() + result[k].Imaginary() * result[k].Imaginary()) }; + EXPECT_NEAR(mag, 1.0f, 1e-3f); + } +} + +TEST_F(TestFastFourierTransformFloat, parseval_energy_preserved) +{ + constexpr std::array signal{ 0.2f, -0.4f, 0.6f, 0.1f, -0.3f, 0.5f, -0.1f, 0.8f }; + timeDomain.resize(Length); + for (std::size_t n{ 0 }; n < Length; ++n) + timeDomain[n] = signal[n]; + + float timePower{ 0.0f }; + for (float s : signal) + timePower += s * s; + + auto& result = fft->Forward(timeDomain); + + float freqPower{ 0.0f }; + for (std::size_t k{ 0 }; k < Length; ++k) + freqPower += result[k].Real() * result[k].Real() + result[k].Imaginary() * result[k].Imaginary(); + freqPower /= static_cast(Length); + + EXPECT_NEAR(freqPower, timePower, 1e-3f); +} + +TEST_F(TestFastFourierTransformFloat, linearity_forward_transform) +{ + constexpr float a{ 3.0f }; + constexpr float b{ -2.0f }; + constexpr std::array x1{ 1.0f, 0.5f, -0.5f, -1.0f, -0.5f, 0.5f, 1.0f, 0.5f }; + constexpr std::array x2{ 0.3f, -0.1f, 0.4f, -0.2f, 0.5f, -0.3f, 0.1f, -0.4f }; + + typename VectorReal::template WithMaxSize in1{}; + typename VectorReal::template WithMaxSize in2{}; + typename VectorReal::template WithMaxSize inCombined{}; + in1.resize(Length); + in2.resize(Length); + inCombined.resize(Length); + + for (std::size_t n{ 0 }; n < Length; ++n) + { + in1[n] = x1[n]; + in2[n] = x2[n]; + inCombined[n] = a * x1[n] + b * x2[n]; + } + + auto& s1{ fft->Forward(in1) }; + typename VectorComplex::template WithMaxSize s1Copy{}; + s1Copy.resize(Length); + for (std::size_t k{ 0 }; k < Length; ++k) + s1Copy[k] = s1[k]; + + auto& s2{ fft->Forward(in2) }; + typename VectorComplex::template WithMaxSize s2Copy{}; + s2Copy.resize(Length); + for (std::size_t k{ 0 }; k < Length; ++k) + s2Copy[k] = s2[k]; + + auto& sCombined{ fft->Forward(inCombined) }; + + for (std::size_t k{ 0 }; k < Length; ++k) + { + EXPECT_NEAR(sCombined[k].Real(), a * s1Copy[k].Real() + b * s2Copy[k].Real(), 1e-3f); + EXPECT_NEAR(sCombined[k].Imaginary(), a * s1Copy[k].Imaginary() + b * s2Copy[k].Imaginary(), 1e-3f); + } +} + +TEST_F(TestFastFourierTransformFloat, inverse_of_known_spectrum_recovers_time_domain) +{ + constexpr std::array expected{ 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f }; + + frequencyDomain.resize(Length); + for (std::size_t k{ 0 }; k < Length; ++k) + frequencyDomain[k] = math::Complex{ 1.0f, 0.0f }; + + auto& result = fft->Inverse(frequencyDomain); + + for (std::size_t n{ 0 }; n < Length; ++n) + EXPECT_NEAR(result[n], expected[n], 1e-3f); +} + +TEST_F(TestFastFourierTransformFloat, all_bins_match_direct_dft) +{ + constexpr std::array signal{ 0.1f, 0.3f, -0.2f, 0.5f, 0.4f, -0.1f, 0.0f, 0.2f }; + timeDomain.resize(Length); + for (std::size_t n{ 0 }; n < Length; ++n) + timeDomain[n] = signal[n]; + + auto& result = fft->Forward(timeDomain); + + for (std::size_t k{ 0 }; k < Length; ++k) + { + float refReal{ 0.0f }; + float refImag{ 0.0f }; + for (std::size_t n{ 0 }; n < Length; ++n) + { + float angle{ -2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(Length) }; + refReal += signal[n] * std::cos(angle); + refImag += signal[n] * std::sin(angle); + } + EXPECT_NEAR(result[k].Real(), refReal, 1e-3f); + EXPECT_NEAR(result[k].Imaginary(), refImag, 1e-3f); + } +} From 0f7a2f0b38969c879c525045a34471edf5e7ea35 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 11:30:44 +0000 Subject: [PATCH 3/9] increase coverage for realFFT --- .../test/TestRealFastFourierTransform.cpp | 144 +++++++++++++++--- 1 file changed, 127 insertions(+), 17 deletions(-) diff --git a/numerical/analysis/test/TestRealFastFourierTransform.cpp b/numerical/analysis/test/TestRealFastFourierTransform.cpp index 4f338307..ed317d83 100644 --- a/numerical/analysis/test/TestRealFastFourierTransform.cpp +++ b/numerical/analysis/test/TestRealFastFourierTransform.cpp @@ -69,29 +69,33 @@ TEST_F(TestRealFastFourierTransform, dc_input_has_only_dc_bin) auto& spectrum{ rfft.Forward(input) }; - EXPECT_NEAR(spectrum[0].Real(), static_cast(N), 1e-2f); + EXPECT_NEAR(spectrum[0].Real(), static_cast(N), math::Tolerance()); + EXPECT_NEAR(spectrum[0].Imaginary(), 0.0f, math::Tolerance()); for (std::size_t k{ 1 }; k <= N / 2; ++k) { float mag{ std::sqrt(spectrum[k].Real() * spectrum[k].Real() + spectrum[k].Imaginary() * spectrum[k].Imaginary()) }; - EXPECT_NEAR(mag, 0.0f, 1e-2f); + EXPECT_NEAR(mag, 0.0f, math::Tolerance()); } } -TEST_F(TestRealFastFourierTransform, single_real_sinusoid_hits_one_bin) +TEST_F(TestRealFastFourierTransform, sinusoid_bin_magnitude_exact) { + static constexpr std::size_t k0{ 3 }; input.resize(N); for (std::size_t n{ 0 }; n < N; ++n) - input[n] = std::cos(2.0f * std::numbers::pi_v * 3.0f * static_cast(n) / static_cast(N)); + input[n] = std::cos(2.0f * std::numbers::pi_v * static_cast(k0) * static_cast(n) / static_cast(N)); auto& spectrum{ rfft.Forward(input) }; - float mag3{ std::sqrt(spectrum[3].Real() * spectrum[3].Real() + spectrum[3].Imaginary() * spectrum[3].Imaginary()) }; + float mag{ std::sqrt(spectrum[k0].Real() * spectrum[k0].Real() + spectrum[k0].Imaginary() * spectrum[k0].Imaginary()) }; + EXPECT_NEAR(mag, static_cast(N) / 2.0f, 1e-2f); + for (std::size_t k{ 0 }; k <= N / 2; ++k) { - if (k == 3) + if (k == k0) continue; - float mag{ std::sqrt(spectrum[k].Real() * spectrum[k].Real() + spectrum[k].Imaginary() * spectrum[k].Imaginary()) }; - EXPECT_LT(mag, mag3); + float magOther{ std::sqrt(spectrum[k].Real() * spectrum[k].Real() + spectrum[k].Imaginary() * spectrum[k].Imaginary()) }; + EXPECT_NEAR(magOther, 0.0f, 1e-2f); } } @@ -120,11 +124,13 @@ TEST_F(TestRealFastFourierTransform, matches_reference_complex_fft) } } -TEST_F(TestRealFastFourierTransform, nyquist_and_dc_are_real) +TEST_F(TestRealFastFourierTransform, dc_and_nyquist_bins_are_real) { + static constexpr std::array signal{ 0.1f, 0.3f, -0.2f, 0.5f, 0.4f, -0.1f, 0.0f, 0.2f, + -0.3f, 0.6f, 0.1f, -0.4f, 0.2f, 0.0f, -0.1f, 0.3f }; input.resize(N); for (std::size_t n{ 0 }; n < N; ++n) - input[n] = std::cos(2.0f * std::numbers::pi_v * static_cast(n) / static_cast(N)); + input[n] = signal[n]; auto& spectrum{ rfft.Forward(input) }; @@ -132,7 +138,82 @@ TEST_F(TestRealFastFourierTransform, nyquist_and_dc_are_real) EXPECT_NEAR(spectrum[N / 2].Imaginary(), 0.0f, math::Tolerance()); } -TEST_F(TestRealFastFourierTransform, inverse_is_left_inverse) +TEST_F(TestRealFastFourierTransform, hermitian_symmetry_all_bins) +{ + static constexpr std::array signal{ 0.1f, 0.3f, -0.2f, 0.5f, 0.4f, -0.1f, 0.0f, 0.2f, + -0.3f, 0.6f, 0.1f, -0.4f, 0.2f, 0.0f, -0.1f, 0.3f }; + input.resize(N); + for (std::size_t n{ 0 }; n < N; ++n) + input[n] = signal[n]; + + auto& spectrum{ rfft.Forward(input) }; + + for (std::size_t k{ 1 }; k < N / 2; ++k) + { + float refReal{ 0.0f }; + float refImagNk{ 0.0f }; + for (std::size_t n{ 0 }; n < N; ++n) + { + float angleNk{ -2.0f * std::numbers::pi_v * static_cast(N - k) * static_cast(n) / static_cast(N) }; + refReal += signal[n] * std::cos(angleNk); + refImagNk += signal[n] * std::sin(angleNk); + } + EXPECT_NEAR(spectrum[k].Real(), refReal, 1e-2f); + EXPECT_NEAR(spectrum[k].Imaginary(), -refImagNk, 1e-2f); + } +} + +TEST_F(TestRealFastFourierTransform, parseval_energy_conserved) +{ + static constexpr std::array signal{ 0.5f, -0.3f, 0.1f, 0.8f, -0.6f, 0.2f, 0.0f, -0.4f, + 0.7f, -0.1f, 0.3f, -0.5f, 0.4f, 0.1f, -0.2f, 0.6f }; + input.resize(N); + for (std::size_t n{ 0 }; n < N; ++n) + input[n] = signal[n]; + + auto& spectrum{ rfft.Forward(input) }; + + float timeEnergy{ 0.0f }; + for (std::size_t n{ 0 }; n < N; ++n) + timeEnergy += signal[n] * signal[n]; + + float spectralEnergy{ spectrum[0].Real() * spectrum[0].Real() + spectrum[0].Imaginary() * spectrum[0].Imaginary() }; + for (std::size_t k{ 1 }; k < N / 2; ++k) + spectralEnergy += 2.0f * (spectrum[k].Real() * spectrum[k].Real() + spectrum[k].Imaginary() * spectrum[k].Imaginary()); + spectralEnergy += spectrum[N / 2].Real() * spectrum[N / 2].Real() + spectrum[N / 2].Imaginary() * spectrum[N / 2].Imaginary(); + spectralEnergy /= static_cast(N); + + EXPECT_NEAR(spectralEnergy, timeEnergy, 1e-2f); +} + +TEST_F(TestRealFastFourierTransform, impulse_has_flat_magnitude) +{ + input.resize(N, 0.0f); + input[0] = 1.0f; + + auto& spectrum{ rfft.Forward(input) }; + + for (std::size_t k{ 0 }; k <= N / 2; ++k) + { + float mag{ std::sqrt(spectrum[k].Real() * spectrum[k].Real() + spectrum[k].Imaginary() * spectrum[k].Imaginary()) }; + EXPECT_NEAR(mag, 1.0f, 1e-2f); + } +} + +TEST_F(TestRealFastFourierTransform, zero_input_gives_zero_spectrum) +{ + input.resize(N, 0.0f); + + auto& spectrum{ rfft.Forward(input) }; + + for (std::size_t k{ 0 }; k <= N / 2; ++k) + { + EXPECT_NEAR(spectrum[k].Real(), 0.0f, math::Tolerance()); + EXPECT_NEAR(spectrum[k].Imaginary(), 0.0f, math::Tolerance()); + } +} + +TEST_F(TestRealFastFourierTransform, inverse_reconstructs_original_signal) { static constexpr std::array signal{ 0.5f, -0.3f, 0.1f, 0.8f, -0.6f, 0.2f, 0.0f, -0.4f, 0.7f, -0.1f, 0.3f, -0.5f, 0.4f, 0.1f, -0.2f, 0.6f }; @@ -193,16 +274,45 @@ TEST_F(TestRealFastFourierTransform, linearity_holds) } } -TEST_F(TestRealFastFourierTransform, impulse_has_flat_magnitude) +TEST_F(TestRealFastFourierTransform, determinism_same_input_same_output) { - input.resize(N, 0.0f); - input[0] = 1.0f; + static constexpr std::array signal{ 0.1f, 0.3f, -0.2f, 0.5f, 0.4f, -0.1f, 0.0f, 0.2f, + -0.3f, 0.6f, 0.1f, -0.4f, 0.2f, 0.0f, -0.1f, 0.3f }; - auto& spectrum{ rfft.Forward(input) }; + typename infra::BoundedVector::template WithMaxSize in1{}; + typename infra::BoundedVector::template WithMaxSize in2{}; + in1.resize(N); + in2.resize(N); + for (std::size_t n{ 0 }; n < N; ++n) + { + in1[n] = signal[n]; + in2[n] = signal[n]; + } + + auto& s1{ rfft.Forward(in1) }; + typename infra::BoundedVector>::template WithMaxSize s1Copy{}; + s1Copy.resize(N / 2 + 1); + for (std::size_t k{ 0 }; k <= N / 2; ++k) + s1Copy[k] = s1[k]; + + auto& s2{ rfft.Forward(in2) }; for (std::size_t k{ 0 }; k <= N / 2; ++k) { - float mag{ std::sqrt(spectrum[k].Real() * spectrum[k].Real() + spectrum[k].Imaginary() * spectrum[k].Imaginary()) }; - EXPECT_NEAR(mag, 1.0f, 1e-2f); + EXPECT_FLOAT_EQ(s2[k].Real(), s1Copy[k].Real()); + EXPECT_FLOAT_EQ(s2[k].Imaginary(), s1Copy[k].Imaginary()); } } + +TEST_F(TestRealFastFourierTransform, inverse_of_impulse_spectrum_recovers_impulse) +{ + input.resize(N, 0.0f); + input[0] = 1.0f; + + auto& spectrum{ rfft.Forward(input) }; + auto& recovered{ rfft.Inverse(spectrum) }; + + EXPECT_NEAR(recovered[0], 1.0f, 1e-2f); + for (std::size_t n{ 1 }; n < N; ++n) + EXPECT_NEAR(recovered[n], 0.0f, 1e-2f); +} From cc7c4d5397af55e1c1a15bda9bf9ec0220b6a0f3 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 11:42:25 +0000 Subject: [PATCH 4/9] fix psd --- numerical/analysis/PowerDensitySpectrum.cpp | 1 + numerical/analysis/PowerDensitySpectrum.hpp | 13 +- .../test/TestPowerDensitySpectrum.cpp | 158 +++++++++++++++++- 3 files changed, 161 insertions(+), 11 deletions(-) diff --git a/numerical/analysis/PowerDensitySpectrum.cpp b/numerical/analysis/PowerDensitySpectrum.cpp index 0f3188ea..18c8060f 100644 --- a/numerical/analysis/PowerDensitySpectrum.cpp +++ b/numerical/analysis/PowerDensitySpectrum.cpp @@ -4,4 +4,5 @@ namespace analysis { template class PowerSpectralDensity, test::TwiddleFactorsStub, 50>; + template class PowerSpectralDensity, test::TwiddleFactorsStub, 0>; } diff --git a/numerical/analysis/PowerDensitySpectrum.hpp b/numerical/analysis/PowerDensitySpectrum.hpp index 0c8b50ad..f7e66908 100644 --- a/numerical/analysis/PowerDensitySpectrum.hpp +++ b/numerical/analysis/PowerDensitySpectrum.hpp @@ -46,10 +46,8 @@ namespace analysis void ResetSegment(); private: - static constexpr QNumberType segmentSizeInverted = QNumberType(1.0f / static_cast(SegmentSize)); static constexpr std::size_t overlapSize = (SegmentSize * Overlap) / 100; static constexpr std::size_t step = SegmentSize - overlapSize; - QNumberType frequencyResolution; windowing::Window& window; QNumberType samplingTimeInSeconds; TwiddleFactor twiddleFactors; @@ -65,7 +63,6 @@ namespace analysis windowing::Window& window, QNumberType samplingTimeInSeconds) : window(window) , samplingTimeInSeconds(samplingTimeInSeconds) - , frequencyResolution(QNumberType(samplingTimeInSeconds * segmentSizeInverted)) {} template @@ -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(SegmentSize)); + + ++segmentCount; } + float normalization = math::ToFloat(samplingTimeInSeconds) / + (math::ToFloat(window.Power(SegmentSize)) * static_cast(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, test::TwiddleFactorsStub, 50>; + extern template class PowerSpectralDensity, test::TwiddleFactorsStub, 0>; #endif } diff --git a/numerical/analysis/test/TestPowerDensitySpectrum.cpp b/numerical/analysis/test/TestPowerDensitySpectrum.cpp index 56d270b1..231d3279 100644 --- a/numerical/analysis/test/TestPowerDensitySpectrum.cpp +++ b/numerical/analysis/test/TestPowerDensitySpectrum.cpp @@ -27,8 +27,31 @@ namespace } }; + template + class TestPowerSpectralDensityZeroOverlap + : public ::testing::Test + { + public: + static constexpr std::size_t length = 512; + static constexpr std::size_t overlap = 0; + + using Fft = analysis::test::FftStub; + using Twiddle = analysis::test::TwiddleFactorsStub; + using PowerDensitySpectrum = analysis::PowerSpectralDensity; + + analysis::test::WindowStub window; + T samplingTime = T(1.0f / 48000.0f); + std::optional powerDensitySpectrum; + + void SetUp() override + { + powerDensitySpectrum.emplace(window, samplingTime); + } + }; + using TestedTypes = ::testing::Types; TYPED_TEST_SUITE(TestPowerSpectralDensity, TestedTypes); + TYPED_TEST_SUITE(TestPowerSpectralDensityZeroOverlap, TestedTypes); } TYPED_TEST(TestPowerSpectralDensity, when_input_smaller_than_fft_size_throws_assertion) @@ -40,18 +63,139 @@ TYPED_TEST(TestPowerSpectralDensity, when_input_smaller_than_fft_size_throws_ass EXPECT_DEATH(this->powerDensitySpectrum->Calculate(input), ""); } -TYPED_TEST(TestPowerSpectralDensity, overlapping_segments_are_properly_averaged) +TYPED_TEST(TestPowerSpectralDensity, output_size_equals_half_segment_size_plus_one) +{ + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize input; + for (std::size_t i = 0; i < this->length; ++i) + input.push_back(TypeParam(0.5f)); + + auto& result = this->powerDensitySpectrum->Calculate(input); + + EXPECT_EQ(result.size(), TestFixture::length / 2 + 1); +} + +TYPED_TEST(TestPowerSpectralDensity, non_dc_bins_are_zero_for_stub_fft) { float tolerance = math::Tolerance(); + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize input; + for (std::size_t i = 0; i < this->length; ++i) + input.push_back(TypeParam(0.5f)); + + auto& result = this->powerDensitySpectrum->Calculate(input); + + for (std::size_t k = 1; k < result.size(); ++k) + EXPECT_NEAR(math::ToFloat(result[k]), 0.0f, tolerance); +} + +TYPED_TEST(TestPowerSpectralDensity, dc_bin_is_nonzero_for_stub_fft) +{ + if constexpr (!std::is_floating_point_v) + GTEST_SKIP(); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize input; + for (std::size_t i = 0; i < this->length; ++i) + input.push_back(TypeParam(0.5f)); + + auto& result = this->powerDensitySpectrum->Calculate(input); + + EXPECT_GT(math::ToFloat(result[0]), 0.0f); +} + +TYPED_TEST(TestPowerSpectralDensity, repeated_calculate_with_same_input_is_deterministic) +{ + float tolerance = math::Tolerance(); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize input; + for (std::size_t i = 0; i < this->length; ++i) + input.push_back(TypeParam(0.5f)); + + auto& result1 = this->powerDensitySpectrum->Calculate(input); + std::array snapshot; + for (std::size_t i = 0; i < result1.size(); ++i) + snapshot[i] = math::ToFloat(result1[i]); + + auto& result2 = this->powerDensitySpectrum->Calculate(input); + + for (std::size_t i = 0; i < result2.size(); ++i) + EXPECT_NEAR(math::ToFloat(result2[i]), snapshot[i], tolerance); +} + +TYPED_TEST(TestPowerSpectralDensity, welch_averaging_produces_segment_count_independent_dc_bin) +{ + if constexpr (!std::is_floating_point_v) + GTEST_SKIP(); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize inputOne; + for (std::size_t i = 0; i < this->length; ++i) + inputOne.push_back(TypeParam(0.5f)); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize<3 * TestFixture::length / 2> inputTwo; + for (std::size_t i = 0; i < 3 * this->length / 2; ++i) + inputTwo.push_back(TypeParam(0.5f)); + + auto& r1 = this->powerDensitySpectrum->Calculate(inputOne); + float dc1 = math::ToFloat(r1[0]); + + auto& r2 = this->powerDensitySpectrum->Calculate(inputTwo); + float dc2 = math::ToFloat(r2[0]); + + ASSERT_GT(dc1, 0.0f); + float ratio = dc2 / dc1; + EXPECT_NEAR(ratio, 1.0f, 1e-3f); +} + +TYPED_TEST(TestPowerSpectralDensity, dc_bin_equals_hand_computed_reference_for_single_segment) +{ + if constexpr (!std::is_floating_point_v) + GTEST_SKIP(); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize input; + for (std::size_t i = 0; i < this->length; ++i) + input.push_back(TypeParam(0.5f)); + + auto& result = this->powerDensitySpectrum->Calculate(input); + + constexpr float samplingTimeF = 1.0f / 48000.0f; + constexpr float segSizeInv = 1.0f / static_cast(TestFixture::length); + constexpr float windowPower = 0.25f; + constexpr float magSq = 0.5f * 0.5f; + constexpr float expected = (magSq * segSizeInv) * (samplingTimeF / windowPower); + + EXPECT_NEAR(math::ToFloat(result[0]), expected, expected * 1e-2f); +} + +TYPED_TEST(TestPowerSpectralDensityZeroOverlap, output_size_equals_half_segment_size_plus_one) +{ typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize<2 * TestFixture::length> input; + for (std::size_t i = 0; i < 2 * this->length; ++i) + input.push_back(TypeParam(0.5f)); + + auto& result = this->powerDensitySpectrum->Calculate(input); + + EXPECT_EQ(result.size(), TestFixture::length / 2 + 1); +} + +TYPED_TEST(TestPowerSpectralDensityZeroOverlap, zero_overlap_step_equals_segment_size) +{ + if constexpr (!std::is_floating_point_v) + GTEST_SKIP(); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize inputOne; + for (std::size_t i = 0; i < this->length; ++i) + inputOne.push_back(TypeParam(0.5f)); + + typename TestFixture::PowerDensitySpectrum::VectorReal::template WithMaxSize<2 * TestFixture::length> inputTwo; + for (std::size_t i = 0; i < 2 * this->length; ++i) + inputTwo.push_back(TypeParam(0.5f)); - for (std::size_t i = 0; i < this->length * 2; ++i) - input.push_back((i % 2) ? TypeParam(0.5f) : TypeParam(-0.5f)); + auto& r1 = this->powerDensitySpectrum->Calculate(inputOne); + float dc1 = math::ToFloat(r1[0]); - auto& spectrum1 = this->powerDensitySpectrum->Calculate(input); - auto& spectrum2 = this->powerDensitySpectrum->Calculate(input); + auto& r2 = this->powerDensitySpectrum->Calculate(inputTwo); + float dc2 = math::ToFloat(r2[0]); - for (std::size_t i = 0; i < spectrum1.size(); ++i) - EXPECT_NEAR(math::ToFloat(spectrum1[i]), math::ToFloat(spectrum2[i]), tolerance); + ASSERT_GT(dc1, 0.0f); + float ratio = dc2 / dc1; + EXPECT_NEAR(ratio, 1.0f, 1e-3f); } From cba409ef5ac773ba0b7177729a08fd9b86f6774e Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 11:51:40 +0000 Subject: [PATCH 5/9] improve coverage for dwt tests --- .../test/TestDiscreteWaveletTransform.cpp | 251 +++++++++++++----- 1 file changed, 186 insertions(+), 65 deletions(-) diff --git a/numerical/analysis/test/TestDiscreteWaveletTransform.cpp b/numerical/analysis/test/TestDiscreteWaveletTransform.cpp index 37c925ef..853b0ec2 100644 --- a/numerical/analysis/test/TestDiscreteWaveletTransform.cpp +++ b/numerical/analysis/test/TestDiscreteWaveletTransform.cpp @@ -1,7 +1,8 @@ #include "numerical/analysis/DiscreteWaveletTransform.hpp" #include "numerical/math/Tolerance.hpp" -#include "gmock/gmock.h" -#include +#include "gtest/gtest.h" +#include +#include namespace { @@ -14,13 +15,14 @@ namespace analysis::WaveletFilters db2{ analysis::MakeDaubechies2() }; analysis::DiscreteWaveletTransform dwtHaar{ haar }; + analysis::DiscreteWaveletTransform dwtDb2{ db2 }; using Signal16 = typename infra::BoundedVector::WithMaxSize; using Signal4 = typename infra::BoundedVector::WithMaxSize<4>; }; } -TEST_F(TestDiscreteWaveletTransform, haar_of_constant_puts_energy_in_approx) +TEST_F(TestDiscreteWaveletTransform, haar_constant_signal_all_details_zero_approx_reference) { Signal16 x; x.resize(N, 1.0f); @@ -33,14 +35,96 @@ TEST_F(TestDiscreteWaveletTransform, haar_of_constant_puts_energy_in_approx) std::size_t offset = dwtHaar.LevelOffset(lvl); std::size_t halfLen = N >> (lvl + 1); for (std::size_t i = 0; i < halfLen; ++i) - EXPECT_NEAR(coeffs[offset + i], 0.0f, 1e-5f); + EXPECT_NEAR(coeffs[offset + i], 0.0f, math::Tolerance()); } - std::size_t approxOffset = N - (N >> 3); - float approxEnergy{ 0.0f }; - for (std::size_t i = approxOffset; i < N; ++i) - approxEnergy += coeffs[i] * coeffs[i]; - EXPECT_GT(approxEnergy, 0.0f); + constexpr float kExpectedApprox = 2.8284271247f; + EXPECT_NEAR(coeffs[14], kExpectedApprox, math::Tolerance()); + EXPECT_NEAR(coeffs[15], kExpectedApprox, math::Tolerance()); +} + +TEST_F(TestDiscreteWaveletTransform, haar_single_level_impulse_reference) +{ + analysis::DiscreteWaveletTransform dwt1{ haar }; + + Signal16 x; + x.push_back(1.0f); + for (std::size_t i = 1; i < N; ++i) + x.push_back(0.0f); + + Signal16 coeffs; + dwt1.Forward(x, coeffs); + + constexpr float kInvSqrt2 = 0.7071067811865476f; + EXPECT_NEAR(coeffs[0], kInvSqrt2, math::Tolerance()); + for (std::size_t i = 1; i < N / 2; ++i) + EXPECT_NEAR(coeffs[i], 0.0f, math::Tolerance()); + EXPECT_NEAR(coeffs[N / 2], kInvSqrt2, math::Tolerance()); + for (std::size_t i = N / 2 + 1; i < N; ++i) + EXPECT_NEAR(coeffs[i], 0.0f, math::Tolerance()); +} + +TEST_F(TestDiscreteWaveletTransform, haar_single_level_qmf_reference) +{ + analysis::DiscreteWaveletTransform dwt1{ haar }; + + Signal4 x; + x.push_back(1.0f); + x.push_back(2.0f); + x.push_back(3.0f); + x.push_back(4.0f); + + Signal4 coeffs; + dwt1.Forward(x, coeffs); + + constexpr float kInvSqrt2 = 0.7071067811865476f; + EXPECT_NEAR(coeffs[0], (1.0f - 2.0f) * kInvSqrt2, 1e-5f); + EXPECT_NEAR(coeffs[1], (3.0f - 4.0f) * kInvSqrt2, 1e-5f); + EXPECT_NEAR(coeffs[2], (1.0f + 2.0f) * kInvSqrt2, 1e-5f); + EXPECT_NEAR(coeffs[3], (3.0f + 4.0f) * kInvSqrt2, 1e-5f); +} + +TEST_F(TestDiscreteWaveletTransform, db2_single_level_qmf_reference) +{ + analysis::DiscreteWaveletTransform dwt1{ db2 }; + + Signal4 x; + x.push_back(1.0f); + x.push_back(2.0f); + x.push_back(3.0f); + x.push_back(4.0f); + + Signal4 coeffs; + dwt1.Forward(x, coeffs); + + EXPECT_NEAR(coeffs[0], 0.0f, math::Tolerance()); + EXPECT_NEAR(coeffs[1], -1.4142135624f, math::Tolerance()); + EXPECT_NEAR(coeffs[2], 2.3107890345f, math::Tolerance()); + EXPECT_NEAR(coeffs[3], 4.7602787773f, math::Tolerance()); +} + +TEST_F(TestDiscreteWaveletTransform, zero_signal_produces_zero_coefficients) +{ + Signal16 x; + x.resize(N, 0.0f); + + Signal16 coeffs; + dwtHaar.Forward(x, coeffs); + + for (std::size_t i = 0; i < N; ++i) + EXPECT_NEAR(coeffs[i], 0.0f, math::Tolerance()); +} + +TEST_F(TestDiscreteWaveletTransform, zero_coefficients_reconstruct_zero) +{ + Signal16 coeffs; + coeffs.resize(N, 0.0f); + + Signal16 xRec; + dwtHaar.Inverse(coeffs, xRec); + + for (std::size_t i = 0; i < N; ++i) + EXPECT_NEAR(xRec[i], 0.0f, math::Tolerance()); } TEST_F(TestDiscreteWaveletTransform, haar_detail_captures_step_edge) @@ -54,15 +138,49 @@ TEST_F(TestDiscreteWaveletTransform, haar_detail_captures_step_edge) Signal16 coeffs; dwtHaar.Forward(x, coeffs); + constexpr float kInvSqrt2 = 0.7071067811865476f; std::size_t offset0 = dwtHaar.LevelOffset(0); - float maxDetail{ 0.0f }; - for (std::size_t i = 0; i < N / 2; ++i) + EXPECT_NEAR(coeffs[offset0 + 2], -kInvSqrt2, math::Tolerance()); +} + +TEST_F(TestDiscreteWaveletTransform, haar_energy_preserved_parseval) +{ + Signal16 x; + for (std::size_t i = 0; i < N; ++i) + x.push_back(static_cast(i) * 0.1f + 0.5f); + + Signal16 coeffs; + dwtHaar.Forward(x, coeffs); + + float energyIn{ 0.0f }; + float energyOut{ 0.0f }; + for (std::size_t i = 0; i < N; ++i) + { + energyIn += x[i] * x[i]; + energyOut += coeffs[i] * coeffs[i]; + } + + EXPECT_NEAR(energyIn, energyOut, math::Tolerance()); +} + +TEST_F(TestDiscreteWaveletTransform, db2_energy_preserved_parseval) +{ + Signal16 x; + for (std::size_t i = 0; i < N; ++i) + x.push_back(static_cast(i) * 0.1f); + + Signal16 coeffs; + dwtDb2.Forward(x, coeffs); + + float energyIn{ 0.0f }; + float energyOut{ 0.0f }; + for (std::size_t i = 0; i < N; ++i) { - float v = std::abs(coeffs[offset0 + i]); - if (v > maxDetail) - maxDetail = v; + energyIn += x[i] * x[i]; + energyOut += coeffs[i] * coeffs[i]; } - EXPECT_GT(maxDetail, 1e-3f); + + EXPECT_NEAR(energyIn, energyOut, math::Tolerance()); } TEST_F(TestDiscreteWaveletTransform, perfect_reconstruction_haar) @@ -78,13 +196,11 @@ TEST_F(TestDiscreteWaveletTransform, perfect_reconstruction_haar) dwtHaar.Inverse(coeffs, xRec); for (std::size_t i = 0; i < N; ++i) - EXPECT_NEAR(xRec[i], x[i], 1e-5f); + EXPECT_NEAR(xRec[i], x[i], math::Tolerance()); } TEST_F(TestDiscreteWaveletTransform, perfect_reconstruction_daubechies) { - analysis::DiscreteWaveletTransform dwtDb2{ db2 }; - Signal16 x; for (std::size_t i = 0; i < N; ++i) x.push_back(static_cast(i % 5) * 0.4f - 1.0f); @@ -96,67 +212,57 @@ TEST_F(TestDiscreteWaveletTransform, perfect_reconstruction_daubechies) dwtDb2.Inverse(coeffs, xRec); for (std::size_t i = 0; i < N; ++i) - EXPECT_NEAR(xRec[i], x[i], 1e-4f); + EXPECT_NEAR(xRec[i], x[i], math::Tolerance()); } -TEST_F(TestDiscreteWaveletTransform, energy_is_preserved_orthogonal) +TEST_F(TestDiscreteWaveletTransform, perfect_reconstruction_sweep_haar) { - Signal16 x; - for (std::size_t i = 0; i < N; ++i) - x.push_back(static_cast(i) * 0.1f + 0.5f); - - Signal16 coeffs; - dwtHaar.Forward(x, coeffs); + std::mt19937 rng{ 9876u }; + std::uniform_real_distribution dist{ -5.0f, 5.0f }; - float energyIn{ 0.0f }; - float energyOut{ 0.0f }; - for (std::size_t i = 0; i < N; ++i) + for (int trial = 0; trial < 50; ++trial) { - energyIn += x[i] * x[i]; - energyOut += coeffs[i] * coeffs[i]; - } + Signal16 x; + for (std::size_t i = 0; i < N; ++i) + x.push_back(dist(rng)); - EXPECT_NEAR(energyIn, energyOut, 1e-3f); + Signal16 coeffs; + dwtHaar.Forward(x, coeffs); + Signal16 xRec; + dwtHaar.Inverse(coeffs, xRec); + + for (std::size_t i = 0; i < N; ++i) + EXPECT_NEAR(xRec[i], x[i], math::Tolerance()); + } } -TEST_F(TestDiscreteWaveletTransform, single_level_matches_manual_qmf) +TEST_F(TestDiscreteWaveletTransform, perfect_reconstruction_sweep_daubechies) { - analysis::DiscreteWaveletTransform dwt1{ haar }; - - Signal4 x; - x.push_back(1.0f); - x.push_back(2.0f); - x.push_back(3.0f); - x.push_back(4.0f); + std::mt19937 rng{ 9876u }; + std::uniform_real_distribution dist{ -5.0f, 5.0f }; - Signal4 coeffs; - dwt1.Forward(x, coeffs); + for (int trial = 0; trial < 50; ++trial) + { + Signal16 x; + for (std::size_t i = 0; i < N; ++i) + x.push_back(dist(rng)); - constexpr float inv_sqrt2 = 0.7071067811865476f; - float expCD0 = (1.0f - 2.0f) * inv_sqrt2; - float expCD1 = (3.0f - 4.0f) * inv_sqrt2; - float expCA0 = (1.0f + 2.0f) * inv_sqrt2; - float expCA1 = (3.0f + 4.0f) * inv_sqrt2; + Signal16 coeffs; + dwtDb2.Forward(x, coeffs); + Signal16 xRec; + dwtDb2.Inverse(coeffs, xRec); - EXPECT_NEAR(coeffs[0], expCD0, 1e-5f); - EXPECT_NEAR(coeffs[1], expCD1, 1e-5f); - EXPECT_NEAR(coeffs[2], expCA0, 1e-5f); - EXPECT_NEAR(coeffs[3], expCA1, 1e-5f); + for (std::size_t i = 0; i < N; ++i) + EXPECT_NEAR(xRec[i], x[i], math::Tolerance()); + } } -TEST_F(TestDiscreteWaveletTransform, multilevel_offsets_are_consistent) +TEST_F(TestDiscreteWaveletTransform, level_offset_matches_analytic_formula) { - std::size_t offset0 = dwtHaar.LevelOffset(0); - std::size_t offset1 = dwtHaar.LevelOffset(1); - std::size_t offset2 = dwtHaar.LevelOffset(2); - - EXPECT_EQ(offset0, 0u); - EXPECT_EQ(offset1, N / 2); - EXPECT_EQ(offset2, N / 2 + N / 4); - - EXPECT_EQ(offset1 - offset0, N / 2); - EXPECT_EQ(offset2 - offset1, N / 4); - EXPECT_EQ(N - offset2, N / 4); + EXPECT_EQ(dwtHaar.LevelOffset(0), 0u); + EXPECT_EQ(dwtHaar.LevelOffset(1), N / 2); + EXPECT_EQ(dwtHaar.LevelOffset(2), N / 2 + N / 4); + EXPECT_EQ(dwtHaar.LevelOffset(3), N / 2 + N / 4 + N / 8); } TEST_F(TestDiscreteWaveletTransform, linearity_holds) @@ -185,5 +291,20 @@ TEST_F(TestDiscreteWaveletTransform, linearity_holds) dwtHaar.Forward(xComb, coeffsComb); for (std::size_t i = 0; i < N; ++i) - EXPECT_NEAR(coeffsComb[i], a * coeffs1[i] + b * coeffs2[i], 1e-4f); + EXPECT_NEAR(coeffsComb[i], a * coeffs1[i] + b * coeffs2[i], math::Tolerance()); +} + +TEST_F(TestDiscreteWaveletTransform, determinism_same_input_same_output) +{ + Signal16 x; + for (std::size_t i = 0; i < N; ++i) + x.push_back(static_cast(i) * 0.7f - 3.0f); + + Signal16 coeffs1; + Signal16 coeffs2; + dwtHaar.Forward(x, coeffs1); + dwtHaar.Forward(x, coeffs2); + + for (std::size_t i = 0; i < N; ++i) + EXPECT_FLOAT_EQ(coeffs1[i], coeffs2[i]); } From c10eb992355bf4a00017ba9adee7da2544292c9f Mon Sep 17 00:00:00 2001 From: gfs Date: Sun, 2 Aug 2026 17:25:25 +0200 Subject: [PATCH 6/9] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- numerical/analysis/test/TestDiscreteCosineTransform.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/numerical/analysis/test/TestDiscreteCosineTransform.cpp b/numerical/analysis/test/TestDiscreteCosineTransform.cpp index a93df39b..15cf794d 100644 --- a/numerical/analysis/test/TestDiscreteCosineTransform.cpp +++ b/numerical/analysis/test/TestDiscreteCosineTransform.cpp @@ -34,8 +34,8 @@ namespace class MockFft : public analysis::FastFourierTransform { public: - MOCK_METHOD(VectorComplex&, Forward, (VectorReal& input), (override)); - MOCK_METHOD(VectorReal&, Inverse, (VectorComplex& input), (override)); + MOCK_METHOD(VectorComplex&, Forward, (VectorReal & input), (override)); + MOCK_METHOD(VectorReal&, Inverse, (VectorComplex & input), (override)); }; class TestDiscreteCosineTransform : public ::testing::Test From 27ba96685ff9b9213a32f31182c750ff2a69da14 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 15:30:49 +0000 Subject: [PATCH 7/9] improve hilbert transform --- .../analysis/test/TestHilbertTransform.cpp | 213 +++++++++++++----- 1 file changed, 161 insertions(+), 52 deletions(-) diff --git a/numerical/analysis/test/TestHilbertTransform.cpp b/numerical/analysis/test/TestHilbertTransform.cpp index 1822ad43..8f1cf0a4 100644 --- a/numerical/analysis/test/TestHilbertTransform.cpp +++ b/numerical/analysis/test/TestHilbertTransform.cpp @@ -56,24 +56,11 @@ TEST_F(TestHilbertTransform, cosine_maps_to_sine_imag) { float expectedReal{ std::cos(2.0f * std::numbers::pi_v * f * static_cast(n) / static_cast(N)) }; float expectedImag{ std::sin(2.0f * std::numbers::pi_v * f * static_cast(n) / static_cast(N)) }; - EXPECT_NEAR(a[n].Real(), expectedReal, 1e-2f); - EXPECT_NEAR(a[n].Imaginary(), expectedImag, 1e-2f); + EXPECT_NEAR(a[n].Real(), expectedReal, 1e-4f); + EXPECT_NEAR(a[n].Imaginary(), expectedImag, 1e-4f); } } -TEST_F(TestHilbertTransform, analytic_real_part_equals_input) -{ - static constexpr float f{ 5.0f }; - input.resize(N); - for (std::size_t n{ 0 }; n < N; ++n) - input[n] = std::cos(2.0f * std::numbers::pi_v * f * static_cast(n) / static_cast(N)); - - auto& a{ hilbert.Analytic(input) }; - - for (std::size_t n{ 5 }; n < N - 5; ++n) - EXPECT_NEAR(a[n].Real(), input[n], 1e-2f); -} - TEST_F(TestHilbertTransform, envelope_of_am_signal_is_modulator) { static constexpr float wc{ 2.0f * std::numbers::pi_v * 12.0f / static_cast(N) }; @@ -95,6 +82,26 @@ TEST_F(TestHilbertTransform, envelope_of_am_signal_is_modulator) } } +TEST_F(TestHilbertTransform, analytic_energy_doubles_real_signal_energy) +{ + static constexpr float f{ 5.0f }; + input.resize(N); + for (std::size_t n{ 0 }; n < N; ++n) + input[n] = std::cos(2.0f * std::numbers::pi_v * f * static_cast(n) / static_cast(N)); + + auto& a{ hilbert.Analytic(input) }; + + float realEnergy{ 0.0f }; + float analyticEnergy{ 0.0f }; + for (std::size_t n{ 0 }; n < N; ++n) + { + realEnergy += input[n] * input[n]; + analyticEnergy += a[n].Real() * a[n].Real() + a[n].Imaginary() * a[n].Imaginary(); + } + + EXPECT_NEAR(analyticEnergy, 2.0f * realEnergy, math::Tolerance() * realEnergy); +} + TEST_F(TestHilbertTransform, instantaneous_frequency_of_tone_is_constant) { static constexpr float f0{ 4.0f }; @@ -111,7 +118,7 @@ TEST_F(TestHilbertTransform, instantaneous_frequency_of_tone_is_constant) float phaseNow{ analysis::AnalyticSignalFft::InstantaneousPhase(a[n]) }; float phasePrev{ analysis::AnalyticSignalFft::InstantaneousPhase(a[n - 1]) }; float freq{ analysis::AnalyticSignalFft::InstantaneousFrequency(phaseNow, phasePrev, ts) }; - EXPECT_NEAR(freq, expectedFreq, 1e-2f); + EXPECT_NEAR(freq, expectedFreq, 1e-4f); } } @@ -129,19 +136,21 @@ TEST_F(TestHilbertTransform, instantaneous_frequency_tracks_chirp) auto& a{ hilbert.Analytic(input) }; static constexpr float ts{ 1.0f }; - float prevFreq{ -1.0f }; + float prevFreq{ 0.0f }; + bool first{ true }; for (std::size_t n{ 6 }; n < N - 5; ++n) { float phaseNow{ analysis::AnalyticSignalFft::InstantaneousPhase(a[n]) }; float phasePrev{ analysis::AnalyticSignalFft::InstantaneousPhase(a[n - 1]) }; float freq{ analysis::AnalyticSignalFft::InstantaneousFrequency(phaseNow, phasePrev, ts) }; - if (prevFreq >= 0.0f) + if (!first) EXPECT_GE(freq, prevFreq - 5e-3f); prevFreq = freq; + first = false; } } -TEST_F(TestHilbertTransform, negative_frequencies_are_zeroed) +TEST_F(TestHilbertTransform, negative_frequencies_suppressed) { static constexpr float f{ 5.0f }; input.resize(N); @@ -150,12 +159,7 @@ TEST_F(TestHilbertTransform, negative_frequencies_are_zeroed) auto& a{ hilbert.Analytic(input) }; - float positiveEnergy{ 0.0f }; - for (std::size_t n{ 0 }; n < N; ++n) - positiveEnergy += a[n].Real() * a[n].Real() + a[n].Imaginary() * a[n].Imaginary(); - - typename infra::BoundedVector>::template WithMaxSize spectrum{}; - spectrum.resize(N); + std::array, N> spectrum{}; for (std::size_t k{ 0 }; k < N; ++k) { float re{ 0.0f }; @@ -172,53 +176,158 @@ TEST_F(TestHilbertTransform, negative_frequencies_are_zeroed) for (std::size_t k{ N / 2 + 1 }; k < N; ++k) { float mag{ std::sqrt(spectrum[k].Real() * spectrum[k].Real() + spectrum[k].Imaginary() * spectrum[k].Imaginary()) }; - EXPECT_NEAR(mag, 0.0f, 5e-2f); + EXPECT_NEAR(mag, 0.0f, math::Tolerance()); } } -TEST_F(TestHilbertTransform, fir_matches_fft_in_passband) +TEST_F(TestHilbertTransform, dc_input_has_zero_hilbert) +{ + input.resize(N, 1.0f); + + auto& a{ hilbert.Analytic(input) }; + + for (std::size_t n{ 0 }; n < N; ++n) + EXPECT_NEAR(a[n].Imaginary(), 0.0f, math::Tolerance()); +} + +TEST_F(TestHilbertTransform, nyquist_tone_has_zero_imaginary) { - static constexpr float f{ 8.0f }; - static constexpr std::size_t delay{ 15 }; input.resize(N); for (std::size_t n{ 0 }; n < N; ++n) - input[n] = std::cos(2.0f * std::numbers::pi_v * f * static_cast(n) / static_cast(N)); + input[n] = (n % 2 == 0) ? 1.0f : -1.0f; auto& a{ hilbert.Analytic(input) }; - typename infra::BoundedVector::template WithMaxSize extended{}; - extended.resize(N + delay, 0.0f); for (std::size_t n{ 0 }; n < N; ++n) - extended[n] = input[n]; + EXPECT_NEAR(a[n].Imaginary(), 0.0f, math::Tolerance()); +} - for (std::size_t n{ 10 }; n < 20; ++n) +TEST_F(TestHilbertTransform, zero_signal_produces_zero_analytic) +{ + input.resize(N, 0.0f); + + auto& a{ hilbert.Analytic(input) }; + + for (std::size_t n{ 0 }; n < N; ++n) { - for (std::size_t s{ 0 }; s < n + delay + 1; ++s) - { - float sample{ s < N ? input[s] : 0.0f }; - fir.Filter(sample); - } + EXPECT_NEAR(a[n].Real(), 0.0f, math::Tolerance()); + EXPECT_NEAR(a[n].Imaginary(), 0.0f, math::Tolerance()); + } +} + +TEST_F(TestHilbertTransform, analytic_is_deterministic) +{ + static constexpr float f{ 7.0f }; + input.resize(N); + for (std::size_t n{ 0 }; n < N; ++n) + input[n] = std::cos(2.0f * std::numbers::pi_v * f * static_cast(n) / static_cast(N)); + + auto& a1{ hilbert.Analytic(input) }; + std::array reals1{}; + std::array imags1{}; + for (std::size_t n{ 0 }; n < N; ++n) + { + reals1[n] = a1[n].Real(); + imags1[n] = a1[n].Imaginary(); } - analysis::HilbertFir firFresh{}; - for (std::size_t s{ 0 }; s < delay; ++s) - firFresh.Filter(input[s]); + TwiddleFactors64 twiddle2{}; + analysis::FastFourierTransformRadix2Impl engine2{ twiddle2 }; + analysis::AnalyticSignalFft hilbert2{ engine2 }; - for (std::size_t n{ delay }; n < N - 5; ++n) + typename infra::BoundedVector::template WithMaxSize input2{}; + input2.resize(N); + for (std::size_t n{ 0 }; n < N; ++n) + input2[n] = input[n]; + + auto& a2{ hilbert2.Analytic(input2) }; + + for (std::size_t n{ 0 }; n < N; ++n) { - auto firOut{ firFresh.Filter(input[n]) }; - float fftEnv{ analysis::AnalyticSignalFft::InstantaneousAmplitude(a[n - delay]) }; - float firEnv{ std::sqrt(firOut.Real() * firOut.Real() + firOut.Imaginary() * firOut.Imaginary()) }; - EXPECT_NEAR(firEnv, fftEnv, 1e-1f); + EXPECT_FLOAT_EQ(a2[n].Real(), reals1[n]); + EXPECT_FLOAT_EQ(a2[n].Imaginary(), imags1[n]); } } -TEST_F(TestHilbertTransform, dc_input_has_zero_hilbert) +TEST_F(TestHilbertTransform, analytic_is_linear) { - input.resize(N, 1.0f); + static constexpr float fa{ 3.0f }; + static constexpr float fb{ 7.0f }; + static constexpr float alpha{ 2.5f }; + static constexpr float beta{ -1.5f }; + + typename infra::BoundedVector::template WithMaxSize xa{}; + typename infra::BoundedVector::template WithMaxSize xb{}; + typename infra::BoundedVector::template WithMaxSize xsum{}; + xa.resize(N); + xb.resize(N); + xsum.resize(N); - auto& a{ hilbert.Analytic(input) }; + for (std::size_t n{ 0 }; n < N; ++n) + { + xa[n] = std::cos(2.0f * std::numbers::pi_v * fa * static_cast(n) / static_cast(N)); + xb[n] = std::cos(2.0f * std::numbers::pi_v * fb * static_cast(n) / static_cast(N)); + xsum[n] = alpha * xa[n] + beta * xb[n]; + } + + TwiddleFactors64 twA{}; + analysis::FastFourierTransformRadix2Impl engA{ twA }; + analysis::AnalyticSignalFft hA{ engA }; - for (std::size_t n{ 1 }; n < N - 1; ++n) - EXPECT_NEAR(a[n].Imaginary(), 0.0f, 1e-2f); + TwiddleFactors64 twB{}; + analysis::FastFourierTransformRadix2Impl engB{ twB }; + analysis::AnalyticSignalFft hB{ engB }; + + auto& aa{ hA.Analytic(xa) }; + std::array combinedReal{}; + std::array combinedImag{}; + for (std::size_t n{ 0 }; n < N; ++n) + { + combinedReal[n] = alpha * aa[n].Real(); + combinedImag[n] = alpha * aa[n].Imaginary(); + } + + auto& ab{ hB.Analytic(xb) }; + for (std::size_t n{ 0 }; n < N; ++n) + { + combinedReal[n] += beta * ab[n].Real(); + combinedImag[n] += beta * ab[n].Imaginary(); + } + + auto& asum{ hilbert.Analytic(xsum) }; + + for (std::size_t n{ 2 }; n < N - 2; ++n) + { + EXPECT_NEAR(asum[n].Real(), combinedReal[n], math::Tolerance()); + EXPECT_NEAR(asum[n].Imaginary(), combinedImag[n], math::Tolerance()); + } +} + +TEST_F(TestHilbertTransform, fir_quadrature_at_passband_frequency) +{ + static constexpr float f{ 4.0f }; + static constexpr std::size_t FirTaps{ 31 }; + static constexpr std::size_t centerTap{ (FirTaps - 1) / 2 }; + static constexpr std::size_t warmup{ FirTaps }; + + analysis::HilbertFir firFilter{}; + + for (std::size_t s{ 0 }; s < warmup; ++s) + { + float val{ std::cos(2.0f * std::numbers::pi_v * f * static_cast(s) / static_cast(N)) }; + firFilter.Filter(val); + } + + float maxErr{ 0.0f }; + for (std::size_t s{ warmup }; s < warmup + 20; ++s) + { + float val{ std::cos(2.0f * std::numbers::pi_v * f * static_cast(s) / static_cast(N)) }; + auto out{ firFilter.Filter(val) }; + std::size_t realSample{ s - centerTap }; + float expectedImag{ std::sin(2.0f * std::numbers::pi_v * f * static_cast(realSample) / static_cast(N)) }; + float err{ std::abs(out.Imaginary() - expectedImag) }; + if (err > maxErr) + maxErr = err; + } + EXPECT_NEAR(maxErr, 0.0f, 5e-2f); } From bb7f01154095f49e7b4f58e045f499ccf0a23a4c Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 17:32:50 +0000 Subject: [PATCH 8/9] add goertzel algorithm --- doc/analysis/GoertzelAlgorithm.md | 22 +- numerical/analysis/GoertzelAlgorithm.hpp | 5 +- .../analysis/test/TestGoertzelAlgorithm.cpp | 201 +++++++++++++++--- 3 files changed, 184 insertions(+), 44 deletions(-) diff --git a/doc/analysis/GoertzelAlgorithm.md b/doc/analysis/GoertzelAlgorithm.md index e30d964e..0aed745c 100644 --- a/doc/analysis/GoertzelAlgorithm.md +++ b/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 diff --git a/numerical/analysis/GoertzelAlgorithm.hpp b/numerical/analysis/GoertzelAlgorithm.hpp index 911c02bf..caf199e4 100644 --- a/numerical/analysis/GoertzelAlgorithm.hpp +++ b/numerical/analysis/GoertzelAlgorithm.hpp @@ -73,8 +73,9 @@ namespace analysis template math::Complex GoertzelAlgorithm::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{ real, imag }; } diff --git a/numerical/analysis/test/TestGoertzelAlgorithm.cpp b/numerical/analysis/test/TestGoertzelAlgorithm.cpp index 56c2aae9..ab072ef5 100644 --- a/numerical/analysis/test/TestGoertzelAlgorithm.cpp +++ b/numerical/analysis/test/TestGoertzelAlgorithm.cpp @@ -1,57 +1,90 @@ #include "numerical/analysis/GoertzelAlgorithm.hpp" #include "numerical/math/Tolerance.hpp" #include "gmock/gmock.h" +#include #include #include -#include namespace { class TestGoertzelAlgorithm : public ::testing::Test { - public: + protected: static constexpr std::size_t N = 16; }; } -TEST_F(TestGoertzelAlgorithm, detects_matching_tone) +TEST_F(TestGoertzelAlgorithm, not_ready_before_n_samples) +{ + analysis::GoertzelAlgorithm goertzel{ std::size_t{ 2 }, N }; + + for (std::size_t n = 0; n < N - 1; ++n) + goertzel.Push(0.0f); + + EXPECT_FALSE(goertzel.Ready()); +} + +TEST_F(TestGoertzelAlgorithm, ready_after_n_samples) { analysis::GoertzelAlgorithm goertzel{ std::size_t{ 2 }, N }; + for (std::size_t n = 0; n < N; ++n) + goertzel.Push(0.0f); + + EXPECT_TRUE(goertzel.Ready()); +} + +TEST_F(TestGoertzelAlgorithm, magnitude_matches_dft_bin_on_tone) +{ + constexpr std::size_t k{ 2 }; + analysis::GoertzelAlgorithm goertzel{ k, N }; + for (std::size_t n = 0; n < N; ++n) { - float sample{ std::cos(2.0f * std::numbers::pi_v * 2.0f * static_cast(n) / static_cast(N)) }; + float sample{ std::cos(2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(N)) }; goertzel.Push(sample); } - EXPECT_TRUE(goertzel.Ready()); EXPECT_NEAR(goertzel.Magnitude(), static_cast(N) / 2.0f, 1e-2f); } -TEST_F(TestGoertzelAlgorithm, rejects_off_bin_tone) +TEST_F(TestGoertzelAlgorithm, magnitude_agrees_with_direct_dft_magnitude) { - analysis::GoertzelAlgorithm goertzel{ std::size_t{ 2 }, N }; + constexpr std::size_t k{ 5 }; + analysis::GoertzelAlgorithm goertzel{ k, N }; + + std::array signal{}; + for (std::size_t n = 0; n < N; ++n) + signal[n] = std::cos(2.0f * std::numbers::pi_v * 2.0f * static_cast(n) / static_cast(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) { - float sample{ std::cos(2.0f * std::numbers::pi_v * 5.0f * static_cast(n) / static_cast(N)) }; - goertzel.Push(sample); + float angle{ -2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(N) }; + refReal += signal[n] * std::cos(angle); + refImag += signal[n] * std::sin(angle); } - EXPECT_LT(goertzel.Magnitude(), 1.0f); + float refMag{ std::sqrt(refReal * refReal + refImag * refImag) }; + EXPECT_NEAR(goertzel.Magnitude(), refMag, 1e-3f); } -TEST_F(TestGoertzelAlgorithm, magnitude_matches_dft_bin) +TEST_F(TestGoertzelAlgorithm, result_matches_direct_dft_bin) { constexpr std::size_t k{ 3 }; analysis::GoertzelAlgorithm goertzel{ k, N }; std::array signal{}; for (std::size_t n = 0; n < N; ++n) - signal[n] = std::cos(2.0f * std::numbers::pi_v * static_cast(n) / static_cast(N)) + 0.5f; + signal[n] = std::cos(2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(N)) + + 0.3f * std::sin(2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(N)); - for (std::size_t n = 0; n < N; ++n) - goertzel.Push(signal[n]); + for (float s : signal) + goertzel.Push(s); float refReal{ 0.0f }; float refImag{ 0.0f }; @@ -67,14 +100,24 @@ TEST_F(TestGoertzelAlgorithm, magnitude_matches_dft_bin) EXPECT_NEAR(result.Imaginary(), refImag, 1e-3f); } -TEST_F(TestGoertzelAlgorithm, coefficient_formula_correct) +TEST_F(TestGoertzelAlgorithm, magnitude_squared_equals_result_norm_squared) { - constexpr std::size_t k{ 3 }; - float expected{ 2.0f * std::cos(2.0f * std::numbers::pi_v * static_cast(k) / static_cast(N)) }; - EXPECT_NEAR(analysis::GoertzelAlgorithm::Coefficient(k, N), expected, math::Tolerance()); + constexpr std::size_t k{ 2 }; + analysis::GoertzelAlgorithm goertzel{ k, N }; + + for (std::size_t n = 0; n < N; ++n) + { + float sample{ std::cos(2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(N)) }; + goertzel.Push(sample); + } + + float mag{ goertzel.Magnitude() }; + math::Complex result{ goertzel.Result() }; + float rSq{ result.Real() * result.Real() + result.Imaginary() * result.Imaginary() }; + EXPECT_NEAR(mag * mag, rSq, 1e-2f); } -TEST_F(TestGoertzelAlgorithm, dc_bin_sums_input) +TEST_F(TestGoertzelAlgorithm, dc_bin_purely_real) { constexpr float c{ 1.5f }; analysis::GoertzelAlgorithm goertzel{ std::size_t{ 0 }, N }; @@ -84,26 +127,80 @@ TEST_F(TestGoertzelAlgorithm, dc_bin_sums_input) math::Complex result{ goertzel.Result() }; EXPECT_NEAR(result.Real(), static_cast(N) * c, 1e-3f); + EXPECT_NEAR(result.Imaginary(), 0.0f, math::Tolerance()); } -TEST_F(TestGoertzelAlgorithm, magnitude_squared_matches_result) +TEST_F(TestGoertzelAlgorithm, nyquist_bin_imaginary_is_zero) { - constexpr std::size_t k{ 2 }; - analysis::GoertzelAlgorithm goertzel{ k, N }; + constexpr std::size_t kNyquist{ N / 2 }; + analysis::GoertzelAlgorithm goertzel{ kNyquist, N }; for (std::size_t n = 0; n < N; ++n) { - float sample{ std::cos(2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(N)) }; + float sample{ std::cos(std::numbers::pi_v * static_cast(n)) }; goertzel.Push(sample); } - float mag{ goertzel.Magnitude() }; math::Complex result{ goertzel.Result() }; - float rSq{ result.Real() * result.Real() + result.Imaginary() * result.Imaginary() }; - EXPECT_NEAR(mag * mag, rSq, 1e-2f); + EXPECT_NEAR(result.Imaginary(), 0.0f, math::Tolerance()); + EXPECT_GT(std::abs(result.Real()), 1.0f); } -TEST_F(TestGoertzelAlgorithm, reset_restarts_block) +TEST_F(TestGoertzelAlgorithm, zero_input_zero_output) +{ + analysis::GoertzelAlgorithm goertzel{ std::size_t{ 3 }, N }; + + for (std::size_t n = 0; n < N; ++n) + goertzel.Push(0.0f); + + EXPECT_NEAR(goertzel.Magnitude(), 0.0f, math::Tolerance()); + math::Complex result{ goertzel.Result() }; + EXPECT_NEAR(result.Real(), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.Imaginary(), 0.0f, math::Tolerance()); +} + +TEST_F(TestGoertzelAlgorithm, off_bin_tone_magnitude_below_half_peak) +{ + constexpr std::size_t kTarget{ 2 }; + constexpr std::size_t kOff{ 5 }; + analysis::GoertzelAlgorithm goertzel{ kTarget, N }; + + for (std::size_t n = 0; n < N; ++n) + { + float sample{ std::cos(2.0f * std::numbers::pi_v * static_cast(kOff) * static_cast(n) / static_cast(N)) }; + goertzel.Push(sample); + } + + EXPECT_LT(goertzel.Magnitude(), static_cast(N) / 4.0f); +} + +TEST_F(TestGoertzelAlgorithm, frequency_hz_constructor_matches_bin_constructor) +{ + constexpr float sampleRate{ 1000.0f }; + constexpr float targetFreq{ 125.0f }; + constexpr std::size_t blockLen{ 8 }; + + analysis::GoertzelAlgorithm byBin{ std::size_t{ 1 }, blockLen }; + analysis::GoertzelAlgorithm byHz{ targetFreq, sampleRate, blockLen }; + + std::array signal{}; + for (std::size_t n = 0; n < blockLen; ++n) + signal[n] = std::cos(2.0f * std::numbers::pi_v * targetFreq * static_cast(n) / sampleRate); + + for (float s : signal) + { + byBin.Push(s); + byHz.Push(s); + } + + EXPECT_NEAR(byBin.Magnitude(), byHz.Magnitude(), math::Tolerance()); + math::Complex rBin{ byBin.Result() }; + math::Complex rHz{ byHz.Result() }; + EXPECT_NEAR(rBin.Real(), rHz.Real(), math::Tolerance()); + EXPECT_NEAR(rBin.Imaginary(), rHz.Imaginary(), math::Tolerance()); +} + +TEST_F(TestGoertzelAlgorithm, reset_restores_not_ready_and_yields_correct_magnitude) { constexpr std::size_t k{ 2 }; analysis::GoertzelAlgorithm goertzel{ k, N }; @@ -123,12 +220,52 @@ TEST_F(TestGoertzelAlgorithm, reset_restarts_block) EXPECT_NEAR(goertzel.Magnitude(), static_cast(N) / 2.0f, 1e-2f); } -TEST_F(TestGoertzelAlgorithm, not_ready_before_n_samples) +TEST_F(TestGoertzelAlgorithm, reset_output_matches_fresh_instance) { - analysis::GoertzelAlgorithm goertzel{ std::size_t{ 2 }, N }; + constexpr std::size_t k{ 3 }; + analysis::GoertzelAlgorithm goertzel{ k, N }; - for (std::size_t n = 0; n < N - 1; ++n) - goertzel.Push(0.0f); + std::array signal{}; + for (std::size_t n = 0; n < N; ++n) + signal[n] = std::cos(2.0f * std::numbers::pi_v * static_cast(k) * static_cast(n) / static_cast(N)); - EXPECT_FALSE(goertzel.Ready()); + for (float s : signal) + goertzel.Push(s); + + goertzel.Reset(); + + for (float s : signal) + goertzel.Push(s); + + analysis::GoertzelAlgorithm fresh{ k, N }; + for (float s : signal) + fresh.Push(s); + + EXPECT_NEAR(goertzel.Magnitude(), fresh.Magnitude(), math::Tolerance()); +} + +TEST_F(TestGoertzelAlgorithm, two_instances_do_not_interfere) +{ + constexpr std::size_t k1{ 2 }; + constexpr std::size_t k2{ 5 }; + analysis::GoertzelAlgorithm g1{ k1, N }; + analysis::GoertzelAlgorithm g2{ k2, N }; + + for (std::size_t n = 0; n < N; ++n) + { + float s1{ std::cos(2.0f * std::numbers::pi_v * static_cast(k1) * static_cast(n) / static_cast(N)) }; + float s2{ std::cos(2.0f * std::numbers::pi_v * static_cast(k2) * static_cast(n) / static_cast(N)) }; + g1.Push(s1); + g2.Push(s2); + } + + EXPECT_NEAR(g1.Magnitude(), static_cast(N) / 2.0f, 1e-2f); + EXPECT_NEAR(g2.Magnitude(), static_cast(N) / 2.0f, 1e-2f); +} + +TEST_F(TestGoertzelAlgorithm, coefficient_formula_correct) +{ + constexpr std::size_t k{ 3 }; + float expected{ 2.0f * std::cos(2.0f * std::numbers::pi_v * static_cast(k) / static_cast(N)) }; + EXPECT_NEAR(analysis::GoertzelAlgorithm::Coefficient(k, N), expected, math::Tolerance()); } From 90ee73966783beac20f78f5d45a0dee2bfb571e5 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 17:36:24 +0000 Subject: [PATCH 9/9] extend tests for decibels and signal detectors --- numerical/analysis/test/TestDecibels.cpp | 7 ++- .../analysis/test/TestSignalDetectors.cpp | 47 +++++++++++++++++++ 2 files changed, 53 insertions(+), 1 deletion(-) diff --git a/numerical/analysis/test/TestDecibels.cpp b/numerical/analysis/test/TestDecibels.cpp index 4c02ec31..f2d0fca6 100644 --- a/numerical/analysis/test/TestDecibels.cpp +++ b/numerical/analysis/test/TestDecibels.cpp @@ -39,6 +39,11 @@ TEST_F(TestDecibels, negative_ratio_returns_floor) EXPECT_NEAR(analysis::ToDecibels(-1.0f), analysis::DecibelFloor::value, math::Tolerance()); } +TEST_F(TestDecibels, tiny_positive_ratio_clamps_to_floor) +{ + EXPECT_NEAR(analysis::ToDecibels(1e-9f), analysis::DecibelFloor::value, math::Tolerance()); +} + TEST_F(TestDecibels, from_decibels_twenty_returns_ten) { EXPECT_NEAR(analysis::FromDecibels(20.0f), 10.0f, math::Tolerance()); @@ -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()); + EXPECT_NEAR(analysis::RippleDb(1.0f, 0.9f), 0.9151f, 1e-3f); } diff --git a/numerical/analysis/test/TestSignalDetectors.cpp b/numerical/analysis/test/TestSignalDetectors.cpp index 5a4565ed..8ddddd9f 100644 --- a/numerical/analysis/test/TestSignalDetectors.cpp +++ b/numerical/analysis/test/TestSignalDetectors.cpp @@ -32,12 +32,31 @@ TEST_F(TestSignalDetectors, peak_hold_decays_when_leaky) EXPECT_NEAR(leaky.Update(0.0f), 0.1f, math::Tolerance()); } +TEST_F(TestSignalDetectors, peak_hold_leaky_new_peak_overtakes_decayed) +{ + analysis::PeakHold leaky{ 0.5f }; + leaky.Update(1.0f); + leaky.Update(0.0f); + EXPECT_NEAR(leaky.Update(0.9f), 0.9f, math::Tolerance()); +} + TEST_F(TestSignalDetectors, peak_hold_tracks_magnitude_not_sign) { EXPECT_NEAR(peak.Update(-0.9f), 0.9f, math::Tolerance()); EXPECT_NEAR(peak.Update(0.1f), 0.9f, math::Tolerance()); } +TEST_F(TestSignalDetectors, peak_hold_reset_determinism) +{ + peak.Update(1.0f); + peak.Update(0.5f); + peak.Reset(); + + analysis::PeakHold fresh{ 1.0f }; + EXPECT_FLOAT_EQ(peak.Update(0.3f), fresh.Update(0.3f)); + EXPECT_FLOAT_EQ(peak.Update(0.7f), fresh.Update(0.7f)); +} + TEST_F(TestSignalDetectors, zero_crossing_counts_sign_changes) { zc.Update(1.0f); @@ -47,6 +66,18 @@ TEST_F(TestSignalDetectors, zero_crossing_counts_sign_changes) EXPECT_EQ(zc.Count(), 3u); } +TEST_F(TestSignalDetectors, zero_crossing_update_returns_true_on_cross) +{ + zc.Update(1.0f); + EXPECT_TRUE(zc.Update(-1.0f)); +} + +TEST_F(TestSignalDetectors, zero_crossing_update_returns_false_when_no_cross) +{ + zc.Update(1.0f); + EXPECT_FALSE(zc.Update(0.5f)); +} + TEST_F(TestSignalDetectors, zero_crossing_hysteresis_rejects_noise) { analysis::ZeroCrossingCounter noisy{ 0.1f }; @@ -56,6 +87,15 @@ TEST_F(TestSignalDetectors, zero_crossing_hysteresis_rejects_noise) EXPECT_EQ(noisy.Count(), 0u); } +TEST_F(TestSignalDetectors, zero_crossing_hysteresis_passes_above_threshold) +{ + analysis::ZeroCrossingCounter filtered{ 0.1f }; + filtered.Update(0.5f); + filtered.Update(-0.5f); + filtered.Update(0.5f); + EXPECT_EQ(filtered.Count(), 2u); +} + TEST_F(TestSignalDetectors, rms_envelope_converges_to_true_rms) { analysis::RmsEnvelope r{ 0.1f }; @@ -81,6 +121,13 @@ TEST_F(TestSignalDetectors, rms_dc_input_returns_magnitude) EXPECT_NEAR(result, 0.4f, 1e-2f); } +TEST_F(TestSignalDetectors, rms_envelope_reset_nonzero_seeds_state) +{ + constexpr float seedValue{ 0.09f }; + rms.Reset(seedValue); + EXPECT_NEAR(rms.Update(std::sqrt(seedValue)), std::sqrt(seedValue), math::Tolerance()); +} + TEST_F(TestSignalDetectors, reset_clears_all_state) { peak.Update(1.0f);