Add fft/spectrum.h: the packed-spectrum view; migrate pvoc, log_mel, test_pvoc (Stage 5)

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 <cassert>
+#include <complex>
+#include <cstddef>
+#include <cstdint>
+#include <type_traits>
+
+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<const float> is a
+    /// read-only view, packed_spectrum<float> 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 <typename Sample>
+    class packed_spectrum {
+      public:
+        using value_type = std::remove_const_t<Sample>;
+        /// 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::is_integral_v<value_type>, std::int64_t, value_type>;
+
+        static_assert(std::is_same_v<value_type, float> || std::is_same_v<value_type, double>
+                          || std::is_same_v<value_type, std::int16_t> || std::is_same_v<value_type, std::int32_t>,
+                      "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 <typename Other>
+            requires(std::is_const_v<Sample> && std::is_same_v<Other, value_type>)
+        constexpr packed_spectrum(const packed_spectrum<Other>& 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<value_type> bin_engineering(std::size_t k) const noexcept
+            requires std::is_floating_point_v<value_type>
+        {
+            assert(k >= 1 && k < m_n / 2);
+            return std::complex<value_type>(m_data[2 * k], -m_data[2 * k + 1]);
+        }
+
+      private:
+        Sample*     m_data;
+        std::size_t m_n;
+    };
+
+} // namespace tap::dsp

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 <vector>
 
 #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<std::ptrdiff_t>(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<const Sample> 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];

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 <algorithm>
 #include <cassert>
 #include <cmath>
+#include <complex>
 #include <cstddef>
 #include <cstdint>
 #include <type_traits>
 #include <vector>
 
 #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<const Sample> analysis(m_frame.data(), m_fft.size());
             const double expected = 2.0 * k_pi * static_cast<double>(m_hop) / static_cast<double>(m_n_size);
             for (int k = 1; k < m_bins - 1; ++k) {
-                const double re    = static_cast<double>(m_frame[static_cast<size_t>(2 * k)]);
-                const double im    = -static_cast<double>(m_frame[static_cast<size_t>(2 * k + 1)]);
-                const double phase = std::atan2(im, re);
+                const std::complex<Sample> bin   = analysis.bin_engineering(static_cast<size_t>(k));
+                const double               re    = static_cast<double>(bin.real());
+                const double               im    = static_cast<double>(bin.imag());
+                const double               phase = std::atan2(im, re);
 
                 double delta = phase - static_cast<double>(m_prev_phase[static_cast<size_t>(k)]) - expected * k;
                 m_prev_phase[static_cast<size_t>(k)] = static_cast<Sample>(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<int>(m_peaks.size());
-            int       lo      = 1;
+            const packed_spectrum<Sample> 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<int>(m_peaks.size());
+            int       lo        = 1;
             for (int pi = 0; pi < n_peaks; ++pi) {
                 const int p         = m_peaks[static_cast<size_t>(pi)];
                 const int hi        = (pi == n_peaks - 1) ? m_bins - 2 : (p + m_peaks[static_cast<size_t>(pi + 1)]) / 2;
@@ -285,8 +292,9 @@ namespace tap::dsp {
                     if (j < 1 || j > m_bins - 2) {
                         continue;
                     }
-                    double re = static_cast<double>(m_frame[static_cast<size_t>(2 * k)]);
-                    double im = -static_cast<double>(m_frame[static_cast<size_t>(2 * k + 1)]);
+                    const std::complex<Sample> bin = analysis.bin_engineering(static_cast<size_t>(k));
+                    double                     re  = static_cast<double>(bin.real());
+                    double                     im  = static_cast<double>(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<size_t>(2 * j)] += static_cast<Sample>(re * cs - im * sn);
-                    m_synth[static_cast<size_t>(2 * j + 1)] -= static_cast<Sample>(re * sn + im * cs); // conjugate back
+                    synthesis.re(static_cast<size_t>(j)) += static_cast<Sample>(re * cs - im * sn);
+                    synthesis.im(static_cast<size_t>(j)) -= static_cast<Sample>(re * sn + im * cs); // conjugate back
                 }
             }
 
@@ -360,12 +368,13 @@ namespace tap::dsp {
                 m_lpc_work[static_cast<size_t>(i)] = static_cast<Sample>(m_lpc_a[static_cast<size_t>(i)]);
             }
             m_fft.forward_inplace(m_lpc_work.data());
-            m_env[0] = 1.0 / std::max(std::abs(static_cast<double>(m_lpc_work[0])), k_env_floor);
+            const packed_spectrum<const Sample> poly(m_lpc_work.data(), m_fft.size());
+            m_env[0] = 1.0 / std::max(std::abs(static_cast<double>(poly.dc())), k_env_floor);
             m_env[static_cast<size_t>(m_bins - 1)] =
-                1.0 / std::max(std::abs(static_cast<double>(m_lpc_work[1])), k_env_floor);
+                1.0 / std::max(std::abs(static_cast<double>(poly.nyquist())), k_env_floor);
             for (int k = 1; k < m_bins - 1; ++k) {
-                const double re               = static_cast<double>(m_lpc_work[static_cast<size_t>(2 * k)]);
-                const double im               = static_cast<double>(m_lpc_work[static_cast<size_t>(2 * k + 1)]);
+                const double re               = static_cast<double>(poly.re(static_cast<size_t>(k)));
+                const double im               = static_cast<double>(poly.im(static_cast<size_t>(k)));
                 m_env[static_cast<size_t>(k)] = 1.0 / std::max(std::sqrt(re * re + im * im), k_env_floor);
             }
         }

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.*")

tests/test_pvoc.cpp

@@ -13,6 +13,7 @@
 
 #include <gtest/gtest.h>
 
+#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<double>(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<const double> spectrum(frame.data(), n);
+        double                                        energy = 0.0;
+        for (size_t k = 1; k < spectrum.num_bins() - 1; ++k) {
             const double f = static_cast<double>(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;