diff --git a/include/tap/dsp/fft/spectrum.h b/include/tap/dsp/fft/spectrum.h new file mode 100644 index 0000000..aa17cdc --- /dev/null +++ b/include/tap/dsp/fft/spectrum.h @@ -0,0 +1,163 @@ +/// @file spectrum.h +/// @brief Non-owning view over a DspTap packed real spectrum. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// One home for the bin arithmetic that every consumer of basic_real_fft used +// to re-derive by hand (data[0] is DC, data[1] is Nyquist, data[2k]/[2k+1] +// are bin k). The view adds no state beyond the pointer and the size, no +// virtuals and no allocation; every accessor is constexpr, noexcept and +// inlines to the same index expression the hand-written code used, so a +// migration onto it is bit-identical by construction. + +#pragma once + +#include +#include +#include +#include +#include + +namespace tap::dsp { + + /// Non-owning view over the packed spectrum that basic_real_fft's forward + /// transform leaves in place and its inverse consumes. + /// + /// The packing, as numbers, for a transform of N real samples (N even; + /// basic_real_fft requires a power of two >= 4). N/2 + 1 bins live in N + /// values a[0..N): + /// - bin[0] = a[0] (DC; its imaginary part is zero and not stored) + /// - bin[N/2] = a[1] (Nyquist; its imaginary part is zero and not stored) + /// - bin[k] = a[2k] + i * a[2k+1] for 1 <= k < N/2 + /// + /// Sign convention: bin[k] = sum_j x[j] * W^(jk) with W = exp(+2*pi*i/N), + /// so the imaginary parts are CONJUGATED relative to the engineering DFT + /// (exp(-2*pi*i/N)). The inverse transform is unnormalized: an in-place + /// round trip needs a 2/N scaling. + /// + /// The native accessors (dc(), nyquist(), re(k), im(k), power(k)) read + /// the spectrum exactly as stored and are the primary interface; spectral + /// products between spectra from this library need no conjugation. + /// bin_engineering(k) is the one convention-flipping accessor, for code + /// that applies textbook phase formulas verbatim; it says so in its name. + /// + /// `Sample` may be const-qualified: packed_spectrum is a + /// read-only view, packed_spectrum also writes through re()/im()/ + /// dc()/nyquist(). A mutable view converts implicitly to the const one. + /// + /// Non-owning: the view holds a pointer and a size and nothing else. It is + /// valid only while the buffer it views is; copying a view (the + /// mutable-to-const conversion included) copies the pointer and the size, + /// never the spectrum. + /// + /// Value types: float and double (the floating profiles) and int16_t and + /// int32_t (the Q15/Q31 fixed-point profiles, which present this same + /// packing). The native accessors are type-agnostic; power() promotes + /// (see power_type); bin_engineering() exists for the floating profiles + /// only, since std::complex over an integer type is unspecified. + /// + /// @pre data points at N values; N is even and >= 2. This is deliberately + /// looser than fft.h's power-of-two >= 4: the packing itself needs + /// only an even N (N = 2 is DC and Nyquist with an empty interior), + /// and the view does not require a power of two. Bin indices are + /// asserted in debug builds and unchecked in release, as everywhere + /// in the Tap libraries. + template + class packed_spectrum { + public: + using value_type = std::remove_const_t; + /// Type of power(): the sample type for the floating profiles (so the + /// product is computed exactly as the hand-written consumers did), and + /// int64_t for the fixed-point ones, where the operands are promoted + /// BEFORE the multiply. Exact for every int16 pair and for every int32 + /// pair except re == im == INT32_MIN, whose sum is 2^63 and overflows; + /// the fixed-point transform's scaling never produces a full-scale + /// pair, so the view does not guard it. + using power_type = std::conditional_t, std::int64_t, value_type>; + + static_assert(std::is_same_v || std::is_same_v + || std::is_same_v || std::is_same_v, + "packed_spectrum views the basic_real_fft profiles: float, double, int16_t (Q15), int32_t (Q31)"); + + /// View N packed values at data. + constexpr packed_spectrum(Sample* data, std::size_t n) noexcept + : m_data(data) + , m_n(n) { + assert(data != nullptr && n >= 2 && (n % 2) == 0); + } + + /// A mutable view converts to the read-only view over the same buffer. + template + requires(std::is_const_v && std::is_same_v) + constexpr packed_spectrum(const packed_spectrum& other) noexcept + : m_data(other.data()) + , m_n(other.size()) {} + + /// Transform size N: the number of packed values. + [[nodiscard]] constexpr std::size_t size() const noexcept { return m_n; } + /// N/2 + 1: the number of distinct bins, DC and Nyquist included. + [[nodiscard]] constexpr std::size_t num_bins() const noexcept { return m_n / 2 + 1; } + /// The packed buffer itself, for handing to a transform. + [[nodiscard]] constexpr Sample* data() const noexcept { return m_data; } + + /// bin[0], the DC term (a[0]); real by construction. + [[nodiscard]] constexpr Sample& dc() const noexcept { return m_data[0]; } + /// bin[N/2], the Nyquist term (a[1]); real by construction. + [[nodiscard]] constexpr Sample& nyquist() const noexcept { return m_data[1]; } + + /// Real part of bin k (a[2k]). + /// @pre 1 <= k < N/2 — DC and Nyquist have their own accessors. + [[nodiscard]] constexpr Sample& re(std::size_t k) const noexcept { + assert(k >= 1 && k < m_n / 2); + return m_data[2 * k]; + } + /// Imaginary part of bin k (a[2k+1]), in the native exp(+i) convention. + /// @pre 1 <= k < N/2 — DC and Nyquist have their own accessors. + [[nodiscard]] constexpr Sample& im(std::size_t k) const noexcept { + assert(k >= 1 && k < m_n / 2); + return m_data[2 * k + 1]; + } + + /// |bin[k]|^2 for any bin, DC and Nyquist included, computed in + /// power_type as re*re + im*im (in that order; consumers that pinned + /// the hand-written product keep their bits). This is |bin[k]|^2 + /// exactly as stored: no factor 2 for the one-sided packing at + /// 1 <= k < N/2 and no 1/N. Over this packing Parseval reads + /// sum_j x[j]^2 = (1/N) * (power(0) + power(N/2) + 2 * sum_{k=1}^{N/2-1} power(k)). + /// @pre 0 <= k <= N/2. + [[nodiscard]] constexpr power_type power(std::size_t k) const noexcept { + assert(k <= m_n / 2); + if (k == 0) { + const power_type dc = m_data[0]; + return dc * dc; + } + if (k == m_n / 2) { + const power_type nyquist = m_data[1]; + return nyquist * nyquist; + } + const power_type re = m_data[2 * k]; + const power_type im = m_data[2 * k + 1]; + return re * re + im * im; + } + + /// Bin k as a complex number in the ENGINEERING convention + /// (exp(-2*pi*i/N)): this CONJUGATES the stored value, returning + /// a[2k] - i * a[2k+1], so textbook phase-vocoder formulas apply + /// verbatim. Not the native reading of the spectrum; a value written + /// back must be conjugated again (see pvoc.h's synthesis pack). + /// Floating profiles only: std::complex over an integer type is + /// unspecified, and the fixed-point consumers work natively. + /// @pre 1 <= k < N/2. + [[nodiscard]] constexpr std::complex bin_engineering(std::size_t k) const noexcept + requires std::is_floating_point_v + { + assert(k >= 1 && k < m_n / 2); + return std::complex(m_data[2 * k], -m_data[2 * k + 1]); + } + + private: + Sample* m_data; + std::size_t m_n; + }; + +} // namespace tap::dsp diff --git a/include/tap/dsp/log_mel.h b/include/tap/dsp/log_mel.h index 70f36fd..92f4a54 100644 --- a/include/tap/dsp/log_mel.h +++ b/include/tap/dsp/log_mel.h @@ -32,9 +32,10 @@ // geometry, never implied. // - FFT: fft_size >= frame, a power of two; the windowed frame occupies // [0, frame) and zeros occupy [frame, fft_size). The transform is the -// unnormalized real DFT of tap::dsp::basic_real_fft (Ooura contract), so -// the power spectrum is |X_k|^2 with X_k = sum x[n] e^{-i 2 pi k n / N} -// up to the sign of the imaginary part, which power discards. +// unnormalized real DFT of tap::dsp::basic_real_fft (the packed spectrum +// defined in fft/spectrum.h), so the power spectrum is |X_k|^2 with +// X_k = sum x[n] e^{-i 2 pi k n / N} up to the sign of the imaginary +// part, which power discards. // - Bin frequencies: f_k = k * sample_rate / fft_size, k in [0, fft_size/2]. // - Mel scale (HTK): mel(f) = 2595 log10(1 + f / 700). Band edges are // bands + 2 points equally spaced in mel between fmin_hz and fmax_hz. @@ -75,6 +76,7 @@ #include #include "tap/dsp/fft.h" +#include "tap/dsp/fft/spectrum.h" namespace tap::dsp { @@ -307,13 +309,15 @@ namespace tap::dsp { } std::fill(m_spec.begin() + static_cast(frame), m_spec.end(), Sample(0)); m_fft.forward_inplace(m_spec.data()); - // Ooura packing: [0] = DC (real), [1] = Nyquist (real), then (re, im) pairs. - m_power[0] = m_spec[0] * m_spec[0]; - m_power[n / 2] = m_spec[1] * m_spec[1]; - for (std::size_t k = 1; k < n / 2; ++k) { - const Sample re = m_spec[2 * k]; - const Sample im = m_spec[2 * k + 1]; - m_power[k] = re * re + im * im; + // The DspTap packed spectrum (see fft/spectrum.h): [0] = DC (real), + // [1] = Nyquist (real), then (re, im) pairs; power() reads all three. + const packed_spectrum spectrum(m_spec.data(), n); + const std::size_t nyquist = spectrum.num_bins() - 1; + + m_power[0] = spectrum.power(0); + m_power[nyquist] = spectrum.power(nyquist); + for (std::size_t k = 1; k < nyquist; ++k) { + m_power[k] = spectrum.power(k); } for (std::size_t b = 0; b < m_g.bands; ++b) { const band& bd = m_bands[b]; diff --git a/include/tap/dsp/pvoc.h b/include/tap/dsp/pvoc.h index 13861de..4253b17 100644 --- a/include/tap/dsp/pvoc.h +++ b/include/tap/dsp/pvoc.h @@ -25,21 +25,24 @@ // Transient smearing on percussive material remains the known trade of the // phase-vocoder class. The transform is tap::dsp::basic_real_fft, so the float // profile rides the vDSP / CMSIS-Helium backends where the build enables them. -// The packed spectrum uses fft.h's conjugated (W = exp(+2*pi*i/N)) convention; -// this class converts to the engineering convention at unpack and back at pack -// so the textbook phase math applies verbatim. +// The spectrum is read through tap::dsp::packed_spectrum (fft/spectrum.h), +// whose native convention is fft.h's conjugated W = exp(+2*pi*i/N); this class +// unpacks through its bin_engineering() accessor (which conjugates) and +// conjugates back at pack so the textbook phase math applies verbatim. #pragma once #include #include #include +#include #include #include #include #include #include "tap/dsp/fft.h" +#include "tap/dsp/fft/spectrum.h" namespace tap::dsp { @@ -214,12 +217,15 @@ namespace tap::dsp { } // per-bin magnitude and instantaneous frequency (engineering-convention - // phases: conjugate fft.h's exp(+i) imaginary parts on unpack) + // phases: bin_engineering() conjugates the packed spectrum's exp(+i) + // imaginary parts on unpack) + const packed_spectrum analysis(m_frame.data(), m_fft.size()); const double expected = 2.0 * k_pi * static_cast(m_hop) / static_cast(m_n_size); for (int k = 1; k < m_bins - 1; ++k) { - const double re = static_cast(m_frame[static_cast(2 * k)]); - const double im = -static_cast(m_frame[static_cast(2 * k + 1)]); - const double phase = std::atan2(im, re); + const std::complex bin = analysis.bin_engineering(static_cast(k)); + const double re = static_cast(bin.real()); + const double im = static_cast(bin.imag()); + const double phase = std::atan2(im, re); double delta = phase - static_cast(m_prev_phase[static_cast(k)]) - expected * k; m_prev_phase[static_cast(k)] = static_cast(phase); @@ -251,10 +257,11 @@ namespace tap::dsp { // synthesis: translate each peak's region rigidly by an integer bin // offset and rotate it by the accumulated residual phase std::fill(m_synth.begin(), m_synth.end(), Sample(0)); - m_synth[0] = m_frame[0]; // DC and Nyquist pass through untouched: they - m_synth[1] = m_frame[1]; // cannot be relocated, and identity stays exact - const int n_peaks = static_cast(m_peaks.size()); - int lo = 1; + const packed_spectrum synthesis(m_synth.data(), m_fft.size()); + synthesis.dc() = analysis.dc(); // DC and Nyquist pass through untouched: they + synthesis.nyquist() = analysis.nyquist(); // cannot be relocated, and identity stays exact + const int n_peaks = static_cast(m_peaks.size()); + int lo = 1; for (int pi = 0; pi < n_peaks; ++pi) { const int p = m_peaks[static_cast(pi)]; const int hi = (pi == n_peaks - 1) ? m_bins - 2 : (p + m_peaks[static_cast(pi + 1)]) / 2; @@ -285,8 +292,9 @@ namespace tap::dsp { if (j < 1 || j > m_bins - 2) { continue; } - double re = static_cast(m_frame[static_cast(2 * k)]); - double im = -static_cast(m_frame[static_cast(2 * k + 1)]); + const std::complex bin = analysis.bin_engineering(static_cast(k)); + double re = static_cast(bin.real()); + double im = static_cast(bin.imag()); if (m_formant) { // keep the envelope in place: excitation from bin k now sits // at bin j, so trade envelope(source) for envelope(target) @@ -295,8 +303,8 @@ namespace tap::dsp { re *= g; im *= g; } - m_synth[static_cast(2 * j)] += static_cast(re * cs - im * sn); - m_synth[static_cast(2 * j + 1)] -= static_cast(re * sn + im * cs); // conjugate back + synthesis.re(static_cast(j)) += static_cast(re * cs - im * sn); + synthesis.im(static_cast(j)) -= static_cast(re * sn + im * cs); // conjugate back } } @@ -360,12 +368,13 @@ namespace tap::dsp { m_lpc_work[static_cast(i)] = static_cast(m_lpc_a[static_cast(i)]); } m_fft.forward_inplace(m_lpc_work.data()); - m_env[0] = 1.0 / std::max(std::abs(static_cast(m_lpc_work[0])), k_env_floor); + const packed_spectrum poly(m_lpc_work.data(), m_fft.size()); + m_env[0] = 1.0 / std::max(std::abs(static_cast(poly.dc())), k_env_floor); m_env[static_cast(m_bins - 1)] = - 1.0 / std::max(std::abs(static_cast(m_lpc_work[1])), k_env_floor); + 1.0 / std::max(std::abs(static_cast(poly.nyquist())), k_env_floor); for (int k = 1; k < m_bins - 1; ++k) { - const double re = static_cast(m_lpc_work[static_cast(2 * k)]); - const double im = static_cast(m_lpc_work[static_cast(2 * k + 1)]); + const double re = static_cast(poly.re(static_cast(k))); + const double im = static_cast(poly.im(static_cast(k))); m_env[static_cast(k)] = 1.0 / std::max(std::sqrt(re * re + im * im), k_env_floor); } } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index eaaa14b..51d2114 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -134,6 +134,7 @@ tap_dsp_add_gtest_executable(tap_dsp_tests test_pvoc.cpp test_quantize.cpp test_sample_traits.cpp + test_spectrum.cpp test_yin.cpp MAIN_FILTER "*-psola_test/1.*:pvoc_test/1.*:yin_test/1.*:log_mel_test/1.*:nn_test/1.*:Kaiser.FastPrototypeMeetsSpec:Kaiser.BalancedPrototypeMeetsSpec:Kaiser.TransparentPrototypeMeetsSpec:Kaiser.EconomyPrototypeMeetsSpec:Kaiser.RationalPhaseCountMeetsSpec:Kaiser.CompensatedSpecsHoldAt16k:MultitoneAnalysis.*") diff --git a/tests/test_pvoc.cpp b/tests/test_pvoc.cpp index 3f42730..a4f7680 100644 --- a/tests/test_pvoc.cpp +++ b/tests/test_pvoc.cpp @@ -13,6 +13,7 @@ #include +#include "tap/dsp/fft/spectrum.h" #include "tap/dsp/pvoc.h" #include "tap/dsp/yin.h" @@ -138,11 +139,12 @@ namespace { frame[i] = w * static_cast(x[x.size() - n + i]); } fft.forward_inplace(frame.data()); - double energy = 0.0; - for (size_t k = 1; k < n / 2; ++k) { + const tap::dsp::packed_spectrum spectrum(frame.data(), n); + double energy = 0.0; + for (size_t k = 1; k < spectrum.num_bins() - 1; ++k) { const double f = static_cast(k) * k_sr / n; if (f >= lo_hz && f <= hi_hz) { - energy += frame[2 * k] * frame[2 * k] + frame[2 * k + 1] * frame[2 * k + 1]; + energy += spectrum.power(k); } } return energy; diff --git a/tests/test_spectrum.cpp b/tests/test_spectrum.cpp new file mode 100644 index 0000000..010b022 --- /dev/null +++ b/tests/test_spectrum.cpp @@ -0,0 +1,327 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Locks down the tap::dsp::packed_spectrum view. Against real transforms +// (float/double): the DC and Nyquist slots, bin k at a[2k] / a[2k+1], the sign +// of im(k) for an on-bin sine under the W = exp(+2*pi*i/N) convention, power() +// at the edges and in the interior, Parseval over the packing, the +// engineering-convention accessor's conjugation, and that a value written +// through the view inverts to the documented tone. Over hand-filled buffers +// (float/double/int16/int32): the slots, write-through, deduction, and +// power()'s promotion before the multiply. + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "tap/dsp/fft.h" +#include "tap/dsp/fft/spectrum.h" + +namespace { + + template + constexpr double k_tolerance = 0.0; + template <> + constexpr double k_tolerance = 1e-12; + template <> + constexpr double k_tolerance = 2e-5; + + template + std::vector forward(std::vector x) { + tap::dsp::basic_real_fft fft(x.size()); + fft.forward_inplace(x.data()); + return x; + } + + template + std::vector tone(std::size_t n, std::size_t bin, bool sine) { + std::vector x(n); + const double w = 2.0 * std::numbers::pi * static_cast(bin) / static_cast(n); + for (std::size_t j = 0; j < n; ++j) { + const double arg = w * static_cast(j); + x[j] = static_cast(sine ? std::sin(arg) : std::cos(arg)); + } + return x; + } + + template + std::vector random_signal(std::size_t n, unsigned seed) { + std::mt19937 gen(seed); + std::uniform_real_distribution dist(-1.0, 1.0); + std::vector x(n); + for (auto& v : x) { + v = static_cast(dist(gen)); + } + return x; + } + + // Transform-backed battery: the two profiles basic_real_fft ships today. + template + class packed_spectrum_test : public ::testing::Test {}; + + using sample_types = ::testing::Types; + TYPED_TEST_SUITE(packed_spectrum_test, sample_types); + + // Hand-filled battery: every value type the view admits, the fixed-point + // profiles included, so Stage 3 does not have to reopen this header. + template + class packed_spectrum_slot_test : public ::testing::Test {}; + + using slot_types = ::testing::Types; + TYPED_TEST_SUITE(packed_spectrum_slot_test, slot_types); + + TYPED_TEST(packed_spectrum_slot_test, SizeAndBinCount) { + std::vector a(1024, TypeParam(0)); + const tap::dsp::packed_spectrum s(a.data(), a.size()); + EXPECT_EQ(s.size(), 1024u); + EXPECT_EQ(s.num_bins(), 513u); + EXPECT_EQ(s.data(), a.data()); + } + + TYPED_TEST(packed_spectrum_test, DcConstantLandsInDc) { + constexpr std::size_t n = 64; + const auto a = forward(std::vector(n, TypeParam(1))); + const double tol = k_tolerance * n; + + const tap::dsp::packed_spectrum s(a.data(), n); + EXPECT_NEAR(s.dc(), static_cast(n), tol); + EXPECT_NEAR(s.nyquist(), 0.0, tol); + for (std::size_t k = 1; k < n / 2; ++k) { + EXPECT_NEAR(s.re(k), 0.0, tol) << "bin " << k; + EXPECT_NEAR(s.im(k), 0.0, tol) << "bin " << k; + } + EXPECT_NEAR(s.power(0), static_cast(n) * n, tol * n); + } + + TYPED_TEST(packed_spectrum_test, AlternatingSignLandsInNyquist) { + constexpr std::size_t n = 64; + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + x[j] = (j % 2 == 0) ? TypeParam(1) : TypeParam(-1); + } + const auto a = forward(x); + const double tol = k_tolerance * n; + + const tap::dsp::packed_spectrum s(a.data(), n); + EXPECT_NEAR(s.dc(), 0.0, tol); + EXPECT_NEAR(s.nyquist(), static_cast(n), tol); + EXPECT_NEAR(s.power(n / 2), static_cast(n) * n, tol * n); + EXPECT_NEAR(s.power(s.num_bins() - 1), static_cast(n) * n, tol * n); + } + + // The documented sign convention (W = exp(+2*pi*i/N)): an on-bin cosine + // is real (+N/2), an on-bin sine is imaginary with a PLUS sign (+N/2) in + // the native reading, and MINUS in the engineering reading. + TYPED_TEST(packed_spectrum_test, OnBinCosineAndSinePinTheSignOfIm) { + constexpr std::size_t n = 128; + constexpr std::size_t bin = 5; + const auto c = forward(tone(n, bin, false)); + const auto s = forward(tone(n, bin, true)); + const double half = static_cast(n) / 2.0; + const double tol = k_tolerance * n; + + const tap::dsp::packed_spectrum cosine(c.data(), n); + const tap::dsp::packed_spectrum sine(s.data(), n); + EXPECT_NEAR(cosine.re(bin), half, tol); + EXPECT_NEAR(cosine.im(bin), 0.0, tol); + EXPECT_NEAR(sine.re(bin), 0.0, tol); + EXPECT_NEAR(sine.im(bin), half, tol); // +N/2, not -N/2 + + const std::complex eng = sine.bin_engineering(bin); + EXPECT_NEAR(eng.real(), 0.0, tol); + EXPECT_NEAR(eng.imag(), -half, tol); // conjugated: the engineering DFT's -N/2 + EXPECT_EQ(eng.real(), sine.re(bin)); + EXPECT_EQ(eng.imag(), -sine.im(bin)); + + // Both tones carry N^2/4 of power in their bin and none elsewhere. + EXPECT_NEAR(cosine.power(bin), half * half, tol * n); + EXPECT_NEAR(sine.power(bin), half * half, tol * n); + for (std::size_t k = 0; k < cosine.num_bins(); ++k) { + if (k == bin) { + continue; + } + EXPECT_NEAR(cosine.power(k), 0.0, tol * n) << "cos leak, bin " << k; + EXPECT_NEAR(sine.power(k), 0.0, tol * n) << "sin leak, bin " << k; + } + } + + TYPED_TEST(packed_spectrum_slot_test, AccessorsAreTheDocumentedSlots) { + // Fill a[i] = i so every slot is distinguishable, and read it back + // through the view: a[0], a[1], a[2k], a[2k+1]. + constexpr std::size_t n = 32; + std::vector a(n); + for (std::size_t i = 0; i < n; ++i) { + a[i] = static_cast(i); + } + const tap::dsp::packed_spectrum s(a.data(), n); + EXPECT_EQ(s.dc(), a[0]); + EXPECT_EQ(s.nyquist(), a[1]); + for (std::size_t k = 1; k < n / 2; ++k) { + EXPECT_EQ(s.re(k), a[2 * k]) << "bin " << k; + EXPECT_EQ(s.im(k), a[2 * k + 1]) << "bin " << k; + EXPECT_EQ(s.power(k), a[2 * k] * a[2 * k] + a[2 * k + 1] * a[2 * k + 1]) << "bin " << k; + } + EXPECT_EQ(s.power(0), a[0] * a[0]); + EXPECT_EQ(s.power(n / 2), a[1] * a[1]); + } + + TYPED_TEST(packed_spectrum_slot_test, MutableViewWritesThrough) { + constexpr std::size_t n = 16; + std::vector a(n, TypeParam(0)); + + const tap::dsp::packed_spectrum s(a.data(), n); + s.dc() = TypeParam(1); + s.nyquist() = TypeParam(2); + s.re(3) = TypeParam(3); + s.im(3) += TypeParam(4); + EXPECT_EQ(a[0], TypeParam(1)); + EXPECT_EQ(a[1], TypeParam(2)); + EXPECT_EQ(a[6], TypeParam(3)); + EXPECT_EQ(a[7], TypeParam(4)); + + // A mutable view converts to the read-only one over the same buffer. + const tap::dsp::packed_spectrum r = s; + EXPECT_EQ(r.data(), a.data()); + EXPECT_EQ(r.re(3), TypeParam(3)); + static_assert(std::is_same_v); + static_assert(std::is_same_v); + } + + TYPED_TEST(packed_spectrum_slot_test, DeductionFollowsThePointer) { + std::vector a(8, TypeParam(0)); + const auto* ca = a.data(); + tap::dsp::packed_spectrum m(a.data(), a.size()); + tap::dsp::packed_spectrum c(ca, a.size()); + static_assert(std::is_same_v>); + static_assert(std::is_same_v>); + EXPECT_EQ(m.num_bins(), 5u); + EXPECT_EQ(c.num_bins(), 5u); + } + + TYPED_TEST(packed_spectrum_slot_test, PowerPromotesBeforeTheMultiply) { + // power_type is the sample type for the floating profiles (the + // hand-written consumers' exact expression) and int64 for the + // fixed-point ones, promoted BEFORE the multiply: an int16 pair of + // 30000 squares to 1.8e9, which int16 arithmetic narrowed to -11776, + // and an int32 pair of 2^30 squares to 2^61, past int32 entirely. + using view = tap::dsp::packed_spectrum; + if constexpr (std::is_integral_v) { + static_assert(std::is_same_v); + } + else { + static_assert(std::is_same_v); + } + + constexpr std::size_t n = 8; + std::vector a(n, TypeParam(0)); + // Chosen with if constexpr, not a ternary: MSVC diagnoses the dead + // TypeParam(1 << 30) branch as a truncating cast under int16 (C4310). + const TypeParam big = [] { + if constexpr (std::is_same_v) { + return TypeParam(1 << 30); + } + else { + return TypeParam(30000); + } + }(); + a[0] = big; + a[1] = big; + a[2] = big; + a[3] = big; + const view s(a.data(), n); + + const std::int64_t b = static_cast(big); + EXPECT_EQ(static_cast(s.power(0)), b * b); + EXPECT_EQ(static_cast(s.power(n / 2)), b * b); + EXPECT_EQ(static_cast(s.power(1)), b * b + b * b); + EXPECT_EQ(s.power(2), typename view::power_type(0)); + } + + TYPED_TEST(packed_spectrum_test, ParsevalHoldsOverThePacking) { + // The power() docstring's formula, with no hidden factor 2 or 1/N: + // sum x^2 = (1/N) (power(0) + power(N/2) + 2 sum_{1..N/2-1} power(k)). + constexpr std::size_t n = 512; + const auto x = random_signal(n, 1234); + + double time_energy = 0.0; + for (std::size_t i = 0; i < n; ++i) { + time_energy += static_cast(x[i]) * static_cast(x[i]); + } + + const auto a = forward(x); + const tap::dsp::packed_spectrum s(a.data(), n); + double freq_energy = static_cast(s.power(0)) + static_cast(s.power(s.num_bins() - 1)); + for (std::size_t k = 1; k < s.num_bins() - 1; ++k) { + freq_energy += 2.0 * static_cast(s.power(k)); + } + freq_energy /= static_cast(n); + + const double tol = std::is_same_v ? 1e-9 : 1e-3; + EXPECT_NEAR(freq_energy, time_energy, tol * time_energy); + } + + TYPED_TEST(packed_spectrum_test, WritesThroughTheViewInvertToTheDocumentedTones) { + // The side pvoc writes: set ONLY im(5) = +N/2 through a mutable view + // and run the real inverse (2/N applied). Under W = exp(+2*pi*i/N) + // that is a sine, not a negative sine; re(5) = N/2 alone is a cosine. + constexpr std::size_t n = 128; + constexpr std::size_t bin = 5; + const double w = 2.0 * std::numbers::pi * static_cast(bin) / static_cast(n); + const double tol = k_tolerance * static_cast(n); + const TypeParam scale = TypeParam(2) / static_cast(n); + + tap::dsp::basic_real_fft fft(n); + + std::vector sine(n, TypeParam(0)); + { + const tap::dsp::packed_spectrum s(sine.data(), n); + s.im(bin) = static_cast(n) / TypeParam(2); + } + fft.inverse_inplace(sine.data()); + + std::vector cosine(n, TypeParam(0)); + { + const tap::dsp::packed_spectrum c(cosine.data(), n); + c.re(bin) = static_cast(n) / TypeParam(2); + } + fft.inverse_inplace(cosine.data()); + + for (std::size_t j = 0; j < n; ++j) { + const double arg = w * static_cast(j); + EXPECT_NEAR(static_cast(sine[j] * scale), std::sin(arg), tol) << "sine, sample " << j; + EXPECT_NEAR(static_cast(cosine[j] * scale), std::cos(arg), tol) << "cosine, sample " << j; + } + } + + TEST(packed_spectrum_contract, AccessorsAreNoexceptAndConstexpr) { + using view = tap::dsp::packed_spectrum; + static_assert(noexcept(std::declval().dc())); + static_assert(noexcept(std::declval().nyquist())); + static_assert(noexcept(std::declval().re(1))); + static_assert(noexcept(std::declval().im(1))); + static_assert(noexcept(std::declval().power(1))); + static_assert(noexcept(std::declval().bin_engineering(1))); + static_assert(noexcept(std::declval().num_bins())); + static_assert(noexcept(std::declval().size())); + + // The whole view is usable in a constant expression: construct it over + // a compile-time buffer and read every accessor back. + constexpr bool reads = [] { + const float a[8] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f}; + const view s(a, 8); + return s.size() == 8 && s.num_bins() == 5 && s.dc() == 1.0f && s.nyquist() == 2.0f && s.re(1) == 3.0f + && s.im(1) == 4.0f && s.power(0) == 1.0f && s.power(2) == 5.0f * 5.0f + 6.0f * 6.0f + && s.power(4) == 4.0f && s.bin_engineering(3).real() == 7.0f && s.bin_engineering(3).imag() == -8.0f; + }(); + static_assert(reads); + SUCCEED(); + } + +} // namespace diff --git a/tools/fingerprint/CMakeLists.txt b/tools/fingerprint/CMakeLists.txt new file mode 100644 index 0000000..e0dfbe5 --- /dev/null +++ b/tools/fingerprint/CMakeLists.txt @@ -0,0 +1,33 @@ +# Same-host A/B fingerprint tool for the DspTap primitives whose output must not +# move when their internals do (pvoc, log_mel). Prints one FNV-1a-64 line per +# case; build it against two trees on ONE host with ONE compiler and flag set, +# and diff the two outputs. The hashes are host + compiler + flags specific +# (pvoc goes through atan2/sin/cos in libm; contraction changes them too), so +# they are never committed as golden values and never asserted in a test. +# Builds on its own: +# +# cmake -B build_fp -S tools/fingerprint -DCMAKE_BUILD_TYPE=Release +# cmake --build build_fp && build_fp/dsptap_fingerprint > after.txt +# +# Copyright 2026 Timothy Place and the DspTap contributors. MIT License. + +cmake_minimum_required(VERSION 3.22) +project(dsptap_fingerprint CXX) + +add_executable(dsptap_fingerprint fingerprint.cpp) + +# pvoc and log_mel pull the shared real FFT (tap::dsp). When this builds as part +# of the DspTap project the target already exists; standalone, pull the repo +# root in directly (its tests stay off: it is not the top-level project). +if (NOT TARGET tap::dsp) + add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/../.. ${CMAKE_CURRENT_BINARY_DIR}/dsptap) +endif () +target_link_libraries(dsptap_fingerprint PRIVATE tap::dsp) +if (TARGET tap_dsp_warnings) + target_link_libraries(dsptap_fingerprint PRIVATE tap_dsp_warnings) +endif () + +set_target_properties(dsptap_fingerprint PROPERTIES + CXX_STANDARD 20 + CXX_STANDARD_REQUIRED ON + CXX_EXTENSIONS OFF) diff --git a/tools/fingerprint/fingerprint.cpp b/tools/fingerprint/fingerprint.cpp new file mode 100644 index 0000000..7936742 --- /dev/null +++ b/tools/fingerprint/fingerprint.cpp @@ -0,0 +1,134 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Same-host A/B fingerprint of the primitives whose output must not move when +// their internals do. Runs tap::dsp::basic_pvoc (float and double, ratio 1.0 +// and 1.5, formant off and on; N = 1024) and tap::dsp::basic_log_mel (float +// and double, plain log and PCEN; the default geometry) over one fixed corpus +// (xorshift32 noise plus two tones, 48000 samples) and prints the FNV-1a-64 +// hash of each case's raw output bytes, one line per case. +// +// HOW TO READ THE NUMBERS. A hash is a function of the source tree AND of the +// host, the compiler, its flags and its libm: pvoc runs atan2/sin/cos in +// double, and floating-point contraction alone changes every pvoc line. So the +// hashes are never golden values, never committed, never asserted in a test +// (a pinned hash would fail on two of the three CI hosts by design). They are +// diffed: build this tool against the tree before a change and the tree after +// it, on the same host with the same compiler and flags, and the twelve lines +// must be identical when the change claims bit identity (the gate for every +// migration of pvoc.h or log_mel.h: the packed-spectrum view, the FFT port's +// routing flip, the engine parameter, the Hann/pi consolidation). +// +// cmake -B build_fp -S tools/fingerprint -DCMAKE_BUILD_TYPE=Release +// cmake --build build_fp && build_fp/dsptap_fingerprint > before.txt +// ... apply the change, rebuild ... +// build_fp/dsptap_fingerprint > after.txt && diff before.txt after.txt +// +// Every case is a chain of noexcept process() calls, so an assert-enabled +// (Debug) build of the same tree exercises the preconditions as well; its +// hashes are comparable only with another Debug build. + +#include +#include +#include +#include +#include + +#include "tap/dsp/log_mel.h" +#include "tap/dsp/pvoc.h" + +namespace { + + constexpr double k_sample_rate = 48000.0; + constexpr std::size_t k_samples = 48000; + constexpr std::size_t k_pvoc_fft = 1024; + + class xorshift32 { + public: + explicit xorshift32(std::uint32_t seed) + : m_s(seed) {} + double next() noexcept { + m_s ^= m_s << 13; + m_s ^= m_s >> 17; + m_s ^= m_s << 5; + return (static_cast(m_s % 65536U) - 32768.0) / 32768.0; + } + + private: + std::uint32_t m_s; + }; + + std::uint64_t fnv1a64(const void* p, std::size_t n) noexcept { + const auto* b = static_cast(p); + std::uint64_t h = 1469598103934665603ULL; + for (std::size_t i = 0; i < n; ++i) { + h ^= b[i]; + h *= 1099511628211ULL; + } + return h; + } + + // The corpus: -20 dB xorshift32 noise under a 220 Hz tone (so pvoc has a + // peak to lock) and an off-bin 1234.5 Hz tone (so it has one to translate + // fractionally), built in double and cast per profile. + std::vector corpus() { + std::vector x(k_samples); + xorshift32 rng(0x2545F491U); + for (std::size_t i = 0; i < k_samples; ++i) { + // Association order is part of the corpus: 2*pi*f*i/sr, left to right. + const double i_d = static_cast(i); + x[i] = 0.1 * rng.next() + 0.5 * std::sin(2.0 * std::numbers::pi * 220.0 * i_d / k_sample_rate) + + 0.3 * std::sin(2.0 * std::numbers::pi * 1234.5 * i_d / k_sample_rate); + } + return x; + } + + void print(const char* primitive, const char* profile, const char* config, std::uint64_t hash) { + std::printf("%-8s %-6s %-22s %016llx\n", primitive, profile, config, static_cast(hash)); + } + + template + void fingerprint_pvoc(const std::vector& x, const char* profile, double ratio, bool formant, + const char* config) { + tap::dsp::basic_pvoc shifter(k_pvoc_fft); + shifter.set_formant(formant); + std::vector out(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + out[i] = shifter.process(static_cast(x[i]), static_cast(ratio)); + } + print("pvoc", profile, config, fnv1a64(out.data(), out.size() * sizeof(Sample))); + } + + template + void fingerprint_log_mel(const std::vector& x, const char* profile, bool pcen, const char* config) { + tap::dsp::log_mel_geometry g; + g.pcen.enabled = pcen; + tap::dsp::basic_log_mel front_end(g); + std::vector in(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + in[i] = static_cast(x[i]); + } + const std::size_t frames = front_end.frames_for(in.size()); + std::vector out(frames * g.bands); + const std::size_t written = front_end.process(in.data(), in.size(), out.data(), frames); + print("log_mel", profile, config, fnv1a64(out.data(), written * g.bands * sizeof(Sample))); + } + + template + void fingerprint_profile(const std::vector& x, const char* profile) { + fingerprint_pvoc(x, profile, 1.0, false, "ratio=1.0 formant=off"); + fingerprint_pvoc(x, profile, 1.5, false, "ratio=1.5 formant=off"); + fingerprint_pvoc(x, profile, 1.0, true, "ratio=1.0 formant=on"); + fingerprint_pvoc(x, profile, 1.5, true, "ratio=1.5 formant=on"); + fingerprint_log_mel(x, profile, false, "log"); + fingerprint_log_mel(x, profile, true, "pcen"); + } + +} // namespace + +int main() { + const std::vector x = corpus(); + fingerprint_profile(x, "float"); + fingerprint_profile(x, "double"); + return 0; +}