diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 69d7992..ae79bfe 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -19,6 +19,8 @@ jobs: runs-on: ${{ matrix.os }} steps: - uses: actions/checkout@v4 + with: + submodules: recursive # the spectral kernels' FFT comes from the dsptap submodule - name: Configure run: cmake -S . -B build -DCMAKE_BUILD_TYPE=Release diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..312d995 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "submodules/dsptap"] + path = submodules/dsptap + url = https://github.com/tap/dsptap diff --git a/CMakeLists.txt b/CMakeLists.txt index 2934059..bacf2f8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,12 +27,21 @@ option(TAPTOOLS_BUILD_TOOLS "Build the offline render tools and the C ABI" ${TAP option(TAPTOOLS_BUILD_BENCH "Build the CPU benchmarks" ${TAPTOOLS_KERNEL_TOP_LEVEL}) option(TAPTOOLS_INSTALL "Generate install and package-export rules" ${TAPTOOLS_KERNEL_TOP_LEVEL}) +# Shared real FFT (DspTap) — the Ooura wrapper used by the spectral kernels +# (conv_engine / stft / nr / spectra), formerly a private radix-2 copy in +# fft.h. Consumed from the tap/dsptap submodule as tap::dsp, which carries the +# split-radix Ooura transform plus the per-target vDSP / CMSIS-Helium float32 +# backends. Pinned as submodules/dsptap (the AmbiTap / MuTap pattern). +add_subdirectory(submodules/dsptap) + add_library(taptools INTERFACE) add_library(TapTools::taptools ALIAS taptools) target_compile_features(taptools INTERFACE cxx_std_20) target_include_directories(taptools INTERFACE $ $) +# tap::dsp brings the real-FFT header and links the Ooura static lib. +target_link_libraries(taptools INTERFACE tap::dsp) if (TAPTOOLS_BUILD_TESTS) enable_testing() @@ -53,10 +62,16 @@ if (TAPTOOLS_INSTALL) include(GNUInstallDirs) install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) - install(TARGETS taptools EXPORT TapToolsKernelTargets) + # The DspTap real FFT is part of taptools' interface (the spectral kernels + # include tap/dsp/fft.h), so its headers and targets join the export set. + install(DIRECTORY submodules/dsptap/include/tap DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) + install(TARGETS taptools tap_dsp tap_dsp_fft EXPORT TapToolsKernelTargets) install(EXPORT TapToolsKernelTargets NAMESPACE TapTools:: DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/TapToolsKernel) + install(FILES submodules/dsptap/third_party/ooura/readme.txt + DESTINATION ${CMAKE_INSTALL_DATADIR}/taptools + RENAME OOURA_LICENSE.txt) configure_package_config_file(cmake/TapToolsKernelConfig.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/TapToolsKernelConfig.cmake diff --git a/bench/CMakeLists.txt b/bench/CMakeLists.txt index 2169bb5..a7c400b 100644 --- a/bench/CMakeLists.txt +++ b/bench/CMakeLists.txt @@ -1,7 +1,7 @@ # CPU benchmarks (kernel-only, no Max dependency). Each bench times its kernel over a config # matrix; the ratchet workflow (compare.py + baselines/) is described in README.md. -foreach (bench svf_bench diode_bench) +foreach (bench svf_bench diode_bench fft_bench) add_executable(${bench} ${bench}.cpp) target_link_libraries(${bench} PRIVATE TapTools::taptools) set_target_properties(${bench} PROPERTIES diff --git a/bench/fft_bench.cpp b/bench/fft_bench.cpp new file mode 100644 index 0000000..512da24 --- /dev/null +++ b/bench/fft_bench.cpp @@ -0,0 +1,176 @@ +/// @file +/// @brief Before/after CPU benchmark for the spectral-kernel FFT swap. +/// @details The spectral kernels (conv_engine / stft / nr / spectra) used to carry a private +/// complex radix-2 Cooley–Tukey FFT; they now share DspTap's real FFT (tap::dsp::real_fft +/// — the split-radix Ooura transform, plus vDSP / CMSIS-Helium float32 backends on Apple +/// and Arm). This times a forward+inverse round trip of each, over the FFT sizes the +/// kernels actually use (a convolution partition of B samples runs a 2B-point transform), +/// so the ratio is the per-transform speedup on this host. +/// +/// The "before" path is a verbatim copy of the retired radix-2 routine (double, full +/// complex, engineering convention, inverse /N) kept HERE only as the benchmark baseline — +/// the kernels no longer contain it. On this Linux/x86 host the win is the split-radix +/// algorithm plus doing a real (not complex) transform; on Apple/Arm the vDSP/CMSIS +/// backends add a further ~3x per transform that this desktop build does not exercise. +/// +/// Usage: fft_bench [--json] [--reps N] +/// @author Timothy Place +// SPDX-License-Identifier: BSD-3-Clause +// Copyright 2026 Timothy Place. + +#include +#include +#include +#include +#include +#include + +#include + +namespace { + + constexpr double k_pi = 3.14159265358979323846; + + // The retired complex radix-2 Cooley–Tukey FFT (baseline only — see file header). + void old_transform(std::vector& re, std::vector& im, bool inverse) { + const int n = static_cast(re.size()); + for (int i = 1, j = 0; i < n; ++i) { + int bit = n >> 1; + for (; j & bit; bit >>= 1) { + j ^= bit; + } + j ^= bit; + if (i < j) { + std::swap(re[i], re[j]); + std::swap(im[i], im[j]); + } + } + for (int len = 2; len <= n; len <<= 1) { + const double ang = 2.0 * k_pi / len * (inverse ? 1.0 : -1.0); + const double wr = std::cos(ang), wi = std::sin(ang); + for (int i = 0; i < n; i += len) { + double cwr = 1.0, cwi = 0.0; + for (int k = 0; k < len / 2; ++k) { + const double vr = re[i + k + len / 2] * cwr - im[i + k + len / 2] * cwi; + const double vi = re[i + k + len / 2] * cwi + im[i + k + len / 2] * cwr; + const double ur = re[i + k]; + const double ui = im[i + k]; + re[i + k] = ur + vr; + im[i + k] = ui + vi; + re[i + k + len / 2] = ur - vr; + im[i + k + len / 2] = ui - vi; + const double ncwr = cwr * wr - cwi * wi; + cwi = cwr * wi + cwi * wr; + cwr = ncwr; + } + } + } + if (inverse) { + for (int i = 0; i < n; ++i) { + re[i] /= n; + im[i] /= n; + } + } + } + + std::vector noise(int n, unsigned seed) { + std::vector v(n); + unsigned s = seed; + for (int i = 0; i < n; ++i) { + s = s * 1664525u + 1013904223u; + v[i] = (static_cast(s) / 2147483648.0 - 1.0); + } + return v; + } + + // Best-of-reps nanoseconds for one forward+inverse round trip at size n. + struct pair_result { + double old_ns; + double new_ns; + double checksum; + }; + + pair_result run_size(int n, int reps, long iters) { + const auto x = noise(n, 4321u + static_cast(n)); + + double best_old = 1e30, best_new = 1e30, checksum = 0.0; + + // Old complex radix-2 round trip. + for (int r = 0; r < reps; ++r) { + std::vector re(n), im(n); + double acc = 0.0; + const auto t0 = std::chrono::steady_clock::now(); + for (long it = 0; it < iters; ++it) { + re.assign(x.begin(), x.end()); + std::fill(im.begin(), im.end(), 0.0); + old_transform(re, im, false); + old_transform(re, im, true); + acc += re[it % n]; + } + const auto t1 = std::chrono::steady_clock::now(); + const double ns = std::chrono::duration(t1 - t0).count() / static_cast(iters); + best_old = std::min(best_old, ns); + checksum += acc; + } + + // DspTap real FFT round trip. + tap::dsp::real_fft fft(static_cast(n)); + for (int r = 0; r < reps; ++r) { + std::vector a(n); + double acc = 0.0; + const auto t0 = std::chrono::steady_clock::now(); + for (long it = 0; it < iters; ++it) { + a.assign(x.begin(), x.end()); + fft.forward_inplace(a.data()); + fft.inverse_inplace(a.data()); + acc += a[it % n]; + } + const auto t1 = std::chrono::steady_clock::now(); + const double ns = std::chrono::duration(t1 - t0).count() / static_cast(iters); + best_new = std::min(best_new, ns); + checksum += acc; + } + + return {best_old, best_new, checksum}; + } + +} // namespace + +int main(int argc, char** argv) { + bool json = false; + int reps = 5; + for (int i = 1; i < argc; ++i) { + if (!std::strcmp(argv[i], "--json")) { + json = true; + } + else if (!std::strcmp(argv[i], "--reps") && i + 1 < argc) { + reps = std::atoi(argv[++i]); + } + } + + const std::vector sizes = {256, 512, 1024, 2048, 4096}; + + if (json) { + std::printf("{\n \"cases\": {\n"); + } + else { + std::printf("%-8s %14s %14s %10s %14s\n", "N", "old ns/xform", "new ns/xform", "speedup", "checksum"); + } + for (size_t s = 0; s < sizes.size(); ++s) { + const int n = sizes[s]; + const long iters = 2'000'000L / n; // constant total work across sizes + const auto res = run_size(n, reps, iters); + if (json) { + std::printf(" \"n%d\": { \"old_ns\": %.1f, \"new_ns\": %.1f, \"speedup\": %.2f }%s\n", n, res.old_ns, + res.new_ns, res.old_ns / res.new_ns, s + 1 < sizes.size() ? "," : ""); + } + else { + std::printf("%-8d %14.1f %14.1f %9.2fx %14.3f\n", n, res.old_ns, res.new_ns, res.old_ns / res.new_ns, + res.checksum); + } + } + if (json) { + std::printf(" }\n}\n"); + } + return 0; +} diff --git a/include/taptools/conv_engine.h b/include/taptools/conv_engine.h index 01616c1..d2c4e09 100644 --- a/include/taptools/conv_engine.h +++ b/include/taptools/conv_engine.h @@ -3,21 +3,25 @@ /// @details The DSP core of tap.convolve~, kept deliberately independent of Max/Min so it can be /// unit-tested directly against a reference convolution (see kernel/tests/conv_engine_test.cpp) /// — buffer~-backed Min objects can't link against the mock test kernel, so the portable -/// engine lives here on its own. Plain C++17: no Jamoma, no min-lib, no external FFT -/// library. +/// engine lives here on its own. Plain C++20: no Jamoma, no min-lib. /// /// Uniformly-partitioned overlap-save: the impulse response is split into equal blocks, -/// each transformed once with an in-house radix-2 FFT; the running output is a -/// frequency-domain multiply-accumulate over a frequency-domain delay line (FDL) of past -/// input spectra. Latency is one partition; everything else is exact linear convolution. +/// each transformed once with the shared DspTap real FFT (tap::dsp::real_fft); the running +/// output is a frequency-domain multiply-accumulate over a frequency-domain delay line +/// (FDL) of past input spectra. Latency is one partition; everything else is exact linear +/// convolution. /// /// True stereo — four IR paths give the full 2×2 response, with the two input channels /// transformed once per block and shared across the paths: /// out_l = in_l ∗ h_LL + in_r ∗ h_RL ; out_r = in_l ∗ h_LR + in_r ∗ h_RR /// Path indexing: path = in_channel * 2 + out_channel, i.e. 0=LL, 1=LR, 2=RL, 3=RR. /// -/// Deferred optimisation: the multiply-accumulate and the IR tables use the full complex -/// spectrum; a real-input half-spectrum (N/2+1 bins) form would halve both CPU and memory. +/// Spectra are stored in the real FFT's packed half-spectrum layout (N reals for an +/// N-point transform: DC in slot 0, Nyquist in slot 1, then re/im pairs for bins +/// 1..N/2-1). That halves both the spectral storage and the multiply-accumulate versus a +/// full-complex form, and — because a product of two Hermitian spectra is Hermitian — the +/// convolution is identical to the full-complex version (the exp(+i) sign convention and +/// the 2/N inverse cancel through the forward→multiply→inverse pipeline). /// @author Timothy Place // SPDX-License-Identifier: BSD-3-Clause // Copyright 2003-2026 Timothy Place. @@ -28,9 +32,10 @@ #include #include #include +#include #include -#include "fft.h" // the shared radix-2 FFT (formerly a private copy of this exact routine) +#include "tap/dsp/fft.h" // the shared DspTap real FFT (split-radix Ooura + vDSP / CMSIS backends) namespace tap::tools { @@ -46,25 +51,22 @@ namespace tap::tools { m_block = std::max(1, blocksize); m_fftsize = 2 * m_block; m_max_parts = std::max(1, max_partitions); + m_fft.emplace(static_cast(m_fftsize)); const int flat = m_max_parts * m_fftsize; for (int path = 0; path < k_paths; ++path) { for (int slot = 0; slot < 2; ++slot) { - m_ir_re[path][slot].assign(flat, 0.0); - m_ir_im[path][slot].assign(flat, 0.0); + m_ir[path][slot].assign(flat, 0.0); } } for (int ch = 0; ch < k_channels; ++ch) { - m_fdl_re[ch].assign(flat, 0.0); - m_fdl_im[ch].assign(flat, 0.0); + m_fdl[ch].assign(flat, 0.0); m_prev[ch].assign(m_block, 0.0); m_inblk[ch].assign(m_block, 0.0); m_outblk[ch].assign(m_block, 0.0); } - m_fre.assign(m_fftsize, 0.0); - m_fim.assign(m_fftsize, 0.0); - m_are.assign(m_fftsize, 0.0); - m_aim.assign(m_fftsize, 0.0); + m_fbuf.assign(m_fftsize, 0.0); + m_acc.assign(m_fftsize, 0.0); m_slot_parts[0] = m_slot_parts[1] = 0; m_active.store(-1, std::memory_order_release); // no IR yet @@ -75,19 +77,18 @@ namespace tap::tools { // Zeroes buffers the audio thread reads; safe to call from a message handler (no reallocation). void clear() { for (int ch = 0; ch < k_channels; ++ch) { - std::fill(m_fdl_re[ch].begin(), m_fdl_re[ch].end(), 0.0); - std::fill(m_fdl_im[ch].begin(), m_fdl_im[ch].end(), 0.0); + std::fill(m_fdl[ch].begin(), m_fdl[ch].end(), 0.0); std::fill(m_prev[ch].begin(), m_prev[ch].end(), 0.0); std::fill(m_inblk[ch].begin(), m_inblk[ch].end(), 0.0); std::fill(m_outblk[ch].begin(), m_outblk[ch].end(), 0.0); } - m_pos = 0; - m_fdl = 0; + m_pos = 0; + m_fdl_pos = 0; } int block_size() const { return m_block; } int max_partitions() const { return m_max_parts; } - bool configured() const { return m_block > 0 && !m_ir_re[0][0].empty(); } + bool configured() const { return m_block > 0 && !m_ir[0][0].empty(); } bool has_ir() const { return m_active.load(std::memory_order_acquire) >= 0; } // Build the four IR paths into the currently-inactive slot and publish it atomically. `paths` holds @@ -106,23 +107,20 @@ namespace tap::tools { for (int path = 0; path < k_paths; ++path) { const float* src = paths[path]; - double* Hre = m_ir_re[path][inactive].data(); - double* Him = m_ir_im[path][inactive].data(); + double* H = m_ir[path][inactive].data(); for (int p = 0; p < P; ++p) { - std::fill(m_fre.begin(), m_fre.end(), 0.0); - std::fill(m_fim.begin(), m_fim.end(), 0.0); + std::fill(m_fbuf.begin(), m_fbuf.end(), 0.0); if (src) { for (int j = 0; j < m_block; ++j) { const int idx = p * m_block + j; if (idx < length) { - m_fre[j] = static_cast(src[idx]) * scale; + m_fbuf[j] = static_cast(src[idx]) * scale; } } } - fft::transform(m_fre, m_fim, false); - std::copy(m_fre.begin(), m_fre.end(), Hre + p * m_fftsize); - std::copy(m_fim.begin(), m_fim.end(), Him + p * m_fftsize); + m_fft->forward_inplace(m_fbuf.data()); + std::copy(m_fbuf.begin(), m_fbuf.end(), H + p * m_fftsize); } } @@ -146,22 +144,35 @@ namespace tap::tools { } private: + // Multiply-accumulate one partition into the packed accumulator: acc += h * x, with DC (slot 0) + // and Nyquist (slot 1) purely real and the interior bins full complex. + void mac_partition(const double* h, const double* x) { + const int half = m_fftsize / 2; + m_acc[0] += h[0] * x[0]; // DC + m_acc[1] += h[1] * x[1]; // Nyquist + for (int k = 1; k < half; ++k) { + const double hr = h[2 * k]; + const double hi = h[2 * k + 1]; + const double xr = x[2 * k]; + const double xi = x[2 * k + 1]; + m_acc[2 * k] += hr * xr - hi * xi; + m_acc[2 * k + 1] += hr * xi + hi * xr; + } + } + // Transform this block's two input frames, store them in the frequency-domain delay line, and form // the four-path multiply-accumulate for both outputs. void process_block() { - const int cur = m_fdl; + const int cur = m_fdl_pos; // 1. Analyse both input channels into the FDL (overlap-save frame = [prev block ; this block]). for (int ch = 0; ch < k_channels; ++ch) { for (int j = 0; j < m_block; ++j) { - m_fre[j] = m_prev[ch][j]; - m_fre[m_block + j] = m_inblk[ch][j]; - m_fim[j] = 0.0; - m_fim[m_block + j] = 0.0; + m_fbuf[j] = m_prev[ch][j]; + m_fbuf[m_block + j] = m_inblk[ch][j]; } - fft::transform(m_fre, m_fim, false); - std::copy(m_fre.begin(), m_fre.end(), m_fdl_re[ch].data() + cur * m_fftsize); - std::copy(m_fim.begin(), m_fim.end(), m_fdl_im[ch].data() + cur * m_fftsize); + m_fft->forward_inplace(m_fbuf.data()); + std::copy(m_fbuf.begin(), m_fbuf.end(), m_fdl[ch].data() + cur * m_fftsize); std::copy(m_inblk[ch].begin(), m_inblk[ch].end(), m_prev[ch].begin()); } @@ -171,64 +182,57 @@ namespace tap::tools { for (int ch = 0; ch < k_channels; ++ch) { std::fill(m_outblk[ch].begin(), m_outblk[ch].end(), 0.0); } - m_fdl = (cur + 1) % m_max_parts; + m_fdl_pos = (cur + 1) % m_max_parts; return; } const int P = m_slot_parts[a]; // 2. For each output channel, accumulate the frequency-domain product over both input channels // and every partition, then inverse-transform and keep the non-aliased second half. + const double scale = 2.0 / static_cast(m_fftsize); // completes the real inverse for (int oc = 0; oc < k_channels; ++oc) { - std::fill(m_are.begin(), m_are.end(), 0.0); - std::fill(m_aim.begin(), m_aim.end(), 0.0); + std::fill(m_acc.begin(), m_acc.end(), 0.0); for (int ic = 0; ic < k_channels; ++ic) { const int path = ic * 2 + oc; - const double* Hre = m_ir_re[path][a].data(); - const double* Him = m_ir_im[path][a].data(); - const double* Xre = m_fdl_re[ic].data(); - const double* Xim = m_fdl_im[ic].data(); + const double* H = m_ir[path][a].data(); + const double* X = m_fdl[ic].data(); for (int p = 0; p < P; ++p) { int slot = cur - p; if (slot < 0) { slot += m_max_parts; } - const double* hr = Hre + p * m_fftsize; - const double* hi = Him + p * m_fftsize; - const double* xr = Xre + slot * m_fftsize; - const double* xi = Xim + slot * m_fftsize; - for (int k = 0; k < m_fftsize; ++k) { - m_are[k] += hr[k] * xr[k] - hi[k] * xi[k]; - m_aim[k] += hr[k] * xi[k] + hi[k] * xr[k]; - } + mac_partition(H + p * m_fftsize, X + slot * m_fftsize); } } - fft::transform(m_are, m_aim, true); + m_fft->inverse_inplace(m_acc.data()); // unscaled; the 2/N below completes the round trip for (int j = 0; j < m_block; ++j) { - m_outblk[oc][j] = m_are[m_block + j]; // overlap-save: discard the aliased first half + m_outblk[oc][j] = m_acc[m_block + j] * scale; // overlap-save: discard the aliased first half } } - m_fdl = (cur + 1) % m_max_parts; + m_fdl_pos = (cur + 1) % m_max_parts; } int m_block{0}; // partition (block) size = latency in samples int m_fftsize{0}; // FFT size = 2 * m_block int m_max_parts{0}; // capacity in partitions + std::optional m_fft; // sized in configure() + std::atomic m_active{-1}; // published IR slot (0/1), or -1 for none. Audio thread reads it. int m_slot_parts[2]{0, 0}; // partition count per slot (kept consistent with m_active) - std::array, 2>, k_paths> m_ir_re, m_ir_im; // IR spectra [path][slot] - std::array, k_channels> m_fdl_re, m_fdl_im; // input FDL [channel] + std::array, 2>, k_paths> m_ir; // packed IR spectra [path][slot] + std::array, k_channels> m_fdl; // packed input FDL [channel] std::array, k_channels> m_prev, m_inblk, m_outblk; - std::vector m_fre, m_fim, m_are, m_aim; // scratch + std::vector m_fbuf, m_acc; // packed scratch (analysis frame / accumulator) - int m_pos{0}; // fill/read index within the current block [0, m_block) - int m_fdl{0}; // ring index of the newest input spectrum + int m_pos{0}; // fill/read index within the current block [0, m_block) + int m_fdl_pos{0}; // ring index of the newest input spectrum }; } // namespace tap::tools diff --git a/include/taptools/fft.h b/include/taptools/fft.h deleted file mode 100644 index dabf7a4..0000000 --- a/include/taptools/fft.h +++ /dev/null @@ -1,71 +0,0 @@ -/// @file -/// @brief Shared in-place radix-2 Cooley–Tukey FFT for the TapTools spectral kernels. -/// @details A single copy of the routine that tap.convolve~ (conv_engine), tap.nr~, and -/// tap.spectra~ all use — it lived, byte-identical, in each of them before being -/// consolidated here. Plain C++17, standard library only: no Jamoma, no external FFT -/// library. The transform is in place; `inverse` divides by N (forward is unscaled), -/// so a forward followed by an inverse reconstructs the input. Size must be a power of -/// two (the callers guarantee this). -/// @author Timothy Place -// SPDX-License-Identifier: BSD-3-Clause -// Copyright 2003-2026 Timothy Place. - -#pragma once - -#include -#include - -namespace tap::tools { - namespace fft { - - inline constexpr double k_pi = 3.14159265358979323846; - - // In-place iterative radix-2 Cooley–Tukey FFT. `inverse` divides by N (forward is unscaled). - // re/im must have the same, power-of-two length. - inline void transform(std::vector& re, std::vector& im, bool inverse) { - const int n = static_cast(re.size()); - - for (int i = 1, j = 0; i < n; ++i) { - int bit = n >> 1; - for (; j & bit; bit >>= 1) { - j ^= bit; - } - j ^= bit; - if (i < j) { - std::swap(re[i], re[j]); - std::swap(im[i], im[j]); - } - } - - for (int len = 2; len <= n; len <<= 1) { - const double ang = 2.0 * k_pi / len * (inverse ? 1.0 : -1.0); - const double wr = std::cos(ang); - const double wi = std::sin(ang); - for (int i = 0; i < n; i += len) { - double cwr = 1.0, cwi = 0.0; - for (int k = 0; k < len / 2; ++k) { - const double vr = re[i + k + len / 2] * cwr - im[i + k + len / 2] * cwi; - const double vi = re[i + k + len / 2] * cwi + im[i + k + len / 2] * cwr; - const double ur = re[i + k]; - const double ui = im[i + k]; - re[i + k] = ur + vr; - im[i + k] = ui + vi; - re[i + k + len / 2] = ur - vr; - im[i + k + len / 2] = ui - vi; - const double ncwr = cwr * wr - cwi * wi; - cwi = cwr * wi + cwi * wr; - cwr = ncwr; - } - } - } - - if (inverse) { - for (int i = 0; i < n; ++i) { - re[i] /= n; - im[i] /= n; - } - } - } - - } // namespace fft -} // namespace tap::tools diff --git a/include/taptools/nr.h b/include/taptools/nr.h index 1aadf35..35de1c3 100644 --- a/include/taptools/nr.h +++ b/include/taptools/nr.h @@ -47,9 +47,11 @@ namespace tap::tools { private: // Attenuate bins whose magnitude is below the threshold (soft-knee downward expansion). + // Operates on the half-spectrum (bins 0..N/2); the real transform mirrors the rest. void gate(std::vector& re, std::vector& im, int N) const { - const double thr = m_threshold; - for (int k = 0; k < N; ++k) { + const double thr = m_threshold; + const int half = N / 2; + for (int k = 0; k <= half; ++k) { const double mag = std::sqrt(re[k] * re[k] + im[k] * im[k]) * (2.0 / N); double gain = 1.0; if (thr > 0.0 && mag < thr) { diff --git a/include/taptools/spectra.h b/include/taptools/spectra.h index 296c5d8..0a834a9 100644 --- a/include/taptools/spectra.h +++ b/include/taptools/spectra.h @@ -46,8 +46,10 @@ namespace tap::tools { } private: - // Reorient the lower-half bins, enforce Hermitian symmetry, and write the result back into - // the scaffold's spectrum buffers for the inverse transform. + // Reorient the half-spectrum bins (0..N/2) and write the result back into the scaffold's + // spectrum buffers for the inverse transform. The real inverse reconstructs the mirrored + // upper half, so no explicit Hermitian mirroring is needed here — only DC and Nyquist must + // stay real. void remap(std::vector& re, std::vector& im, int N) { const int half = N / 2; @@ -64,12 +66,8 @@ namespace tap::tools { } m_oim[0] = 0.0; // DC is real m_oim[half] = 0.0; // Nyquist is real - for (int k = 1; k < half; ++k) { - m_ore[N - k] = m_ore[k]; - m_oim[N - k] = -m_oim[k]; - } - for (int k = 0; k < N; ++k) { + for (int k = 0; k <= half; ++k) { re[k] = m_ore[k]; im[k] = m_oim[k]; } diff --git a/include/taptools/stft.h b/include/taptools/stft.h index d0b0ee6..25748d3 100644 --- a/include/taptools/stft.h +++ b/include/taptools/stft.h @@ -5,9 +5,14 @@ /// buffers, COLA normalisation for perfect reconstruction, and the per-sample pump that /// fires an FFT frame every hop. The only thing that differs between objects is what is /// done to the spectrum between the forward and inverse transforms, so `process` takes -/// that step as a callable: `op(re, im, N)` operates in place on the full N-point -/// complex spectrum. With an identity op, the output reconstructs the input delayed by -/// one FFT frame (the latency). Plain C++17, standard library only. +/// that step as a callable: `op(re, im, N)` operates in place on the half-spectrum — bins +/// 0..N/2, the N/2+1 unique bins of the real transform (im[0] and im[N/2] are zero, DC +/// and Nyquist being real). With an identity op, the output reconstructs the input delayed +/// by one FFT frame (the latency). +/// +/// The transform is the shared DspTap real FFT (tap::dsp::real_fft): half the work of the +/// old complex radix-2 for a real signal, and on Apple / Arm targets it dispatches to the +/// vDSP / CMSIS-Helium backends. Plain C++20, standard library only. /// @author Timothy Place // SPDX-License-Identifier: BSD-3-Clause // Copyright 2003-2026 Timothy Place. @@ -15,16 +20,18 @@ #pragma once #include +#include #include #include -#include "fft.h" +#include "tap/dsp/fft.h" namespace tap::tools { class stft { public: - static constexpr int k_default_overlap = 4; + static constexpr int k_default_overlap = 4; + static constexpr double k_pi = 3.14159265358979323846; // Allocate the window and buffers for the given FFT size (a power of two) and overlap factor. // Resets all running state. Call before processing and whenever the size changes. @@ -35,7 +42,7 @@ namespace tap::tools { m_window.assign(m_fftsize, 0.0); for (int k = 0; k < m_fftsize; ++k) { - m_window[k] = 0.5 - 0.5 * std::cos(2.0 * fft::k_pi * k / m_fftsize); + m_window[k] = 0.5 - 0.5 * std::cos(2.0 * k_pi * k / m_fftsize); } // COLA normalisation: overlap-add `overlap` copies of window^2 and read the steady-state @@ -49,6 +56,8 @@ namespace tap::tools { const double c = cola[m_fftsize / 2]; m_norm = (c > 0.0) ? (1.0 / c) : 1.0; + m_fft.emplace(static_cast(m_fftsize)); + m_spec.assign(m_fftsize, 0.0); m_re.assign(m_fftsize, 0.0); m_im.assign(m_fftsize, 0.0); reset(); @@ -68,7 +77,7 @@ namespace tap::tools { int latency() const { return m_fftsize; } // reconstruction delay with an identity op // Pump n samples through the STFT. `op` is invoked once per hop as op(re, im, N) and may modify - // the spectrum in place; input and output must not alias. + // the half-spectrum (bins 0..N/2) in place; input and output must not alias. template void process(const double* in, double* out, long n, SpectralOp&& op) { for (long i = 0; i < n; ++i) { @@ -88,21 +97,42 @@ namespace tap::tools { private: template void process_frame(SpectralOp&& op) { - const int N = m_fftsize; + const int N = m_fftsize; + const int half = N / 2; - // Assemble the most recent N samples (oldest first) and apply the analysis window. + // Assemble the most recent N samples (oldest first), apply the analysis window, and + // forward-transform in place into DspTap's packed spectrum. for (int k = 0; k < N; ++k) { - m_re[k] = m_inbuf[(m_pos + k) % N] * m_window[k]; - m_im[k] = 0.0; + m_spec[k] = m_inbuf[(m_pos + k) % N] * m_window[k]; } + m_fft->forward_inplace(m_spec.data()); + + // Unpack the packed spectrum into the N/2+1 unique bins for the operation. + m_re[0] = m_spec[0]; // DC (real) + m_im[0] = 0.0; + m_re[half] = m_spec[1]; // Nyquist (real) + m_im[half] = 0.0; + for (int k = 1; k < half; ++k) { + m_re[k] = m_spec[2 * k]; + m_im[k] = m_spec[2 * k + 1]; + } + + op(m_re, m_im, N); // in-place spectral operation on the half-spectrum (bins 0..N/2) - fft::transform(m_re, m_im, false); - op(m_re, m_im, N); // in-place spectral operation on the full N-point spectrum - fft::transform(m_re, m_im, true); + // Repack for the inverse transform. + m_spec[0] = m_re[0]; + m_spec[1] = m_re[half]; + for (int k = 1; k < half; ++k) { + m_spec[2 * k] = m_re[k]; + m_spec[2 * k + 1] = m_im[k]; + } + m_fft->inverse_inplace(m_spec.data()); // unscaled — the 2/N is folded into the overlap-add below - // Synthesis window + overlap-add into the output accumulator. + // Synthesis window + overlap-add into the output accumulator. The real inverse is + // unnormalised, so the 2/N that completes the round trip rides in here with the COLA norm. + const double scale = 2.0 / static_cast(N); for (int k = 0; k < N; ++k) { - m_outbuf[(m_pos + k) % N] += m_re[k] * m_window[k] * m_norm; + m_outbuf[(m_pos + k) % N] += m_spec[k] * scale * m_window[k] * m_norm; } } @@ -114,11 +144,13 @@ namespace tap::tools { int m_pos{0}; int m_hopcount{0}; - std::vector m_window; // Hann (analysis and synthesis) - std::vector m_inbuf; // circular input buffer (length N) - std::vector m_outbuf; // circular overlap-add accumulator (length N) - std::vector m_re; // scratch FFT real - std::vector m_im; // scratch FFT imag + std::optional m_fft; // sized in configure() + std::vector m_window; // Hann (analysis and synthesis) + std::vector m_inbuf; // circular input buffer (length N) + std::vector m_outbuf; // circular overlap-add accumulator (length N) + std::vector m_spec; // packed FFT scratch (time in / spectrum / time out) + std::vector m_re; // unpacked half-spectrum real (bins 0..N/2) + std::vector m_im; // unpacked half-spectrum imag (bins 0..N/2) }; } // namespace tap::tools diff --git a/include/taptools/taptools.h b/include/taptools/taptools.h index 7ef5f30..5086fea 100644 --- a/include/taptools/taptools.h +++ b/include/taptools/taptools.h @@ -11,7 +11,6 @@ #include "bridged_t.h" #include "conv_engine.h" #include "diode_ladder.h" -#include "fft.h" #include "grm_comb.h" #include "grm_pitchaccum.h" #include "ladder.h" diff --git a/submodules/dsptap b/submodules/dsptap new file mode 160000 index 0000000..9cbfbca --- /dev/null +++ b/submodules/dsptap @@ -0,0 +1 @@ +Subproject commit 9cbfbca6670bdc929bab2804c713b60dad8b1c0a diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 452701c..ea97a18 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -14,7 +14,6 @@ FetchContent_MakeAvailable(Catch2) add_executable(taptools_kernel_tests autowah_test.cpp diode_ladder_test.cpp - fft_test.cpp grm_comb_test.cpp nr_test.cpp spectra_test.cpp diff --git a/tests/fft_test.cpp b/tests/fft_test.cpp deleted file mode 100644 index 5ea7a30..0000000 --- a/tests/fft_test.cpp +++ /dev/null @@ -1,122 +0,0 @@ -/// @file -/// @brief Unit tests for the shared radix-2 FFT (tap::tools::fft::transform). -/// @details Checks the forward transform against a naive O(N^2) DFT, forward→inverse -/// reconstruction, linearity of the convention (forward unscaled, inverse /N), and a -/// couple of closed-form cases (impulse → flat spectrum, single cosine → two bins). -// SPDX-License-Identifier: BSD-3-Clause -// Copyright 2003-2026 Timothy Place. - -#include -#include - -#include -#include - -namespace { - - constexpr double k_pi = 3.14159265358979323846; - - // Naive O(N^2) DFT reference, same sign convention as tap::tools::fft (forward: exp(-i...)). - void naive_dft(const std::vector& re, const std::vector& im, std::vector& ore, - std::vector& oim) { - const int n = static_cast(re.size()); - ore.assign(n, 0.0); - oim.assign(n, 0.0); - for (int k = 0; k < n; ++k) { - double sr = 0.0, si = 0.0; - for (int t = 0; t < n; ++t) { - const double ang = -2.0 * k_pi * k * t / n; - const double c = std::cos(ang), s = std::sin(ang); - sr += re[t] * c - im[t] * s; - si += re[t] * s + im[t] * c; - } - ore[k] = sr; - oim[k] = si; - } - } - - std::vector noise(int n, unsigned seed) { - std::vector v(n); - unsigned s = seed; - for (int i = 0; i < n; ++i) { - s = s * 1664525u + 1013904223u; - v[i] = (static_cast(s) / 4294967295.0) * 2.0 - 1.0; - } - return v; - } - - double max_abs_diff(const std::vector& a, const std::vector& b) { - double m = 0.0; - for (size_t i = 0; i < a.size(); ++i) { - m = std::max(m, std::abs(a[i] - b[i])); - } - return m; - } - -} // namespace - -SCENARIO("the forward FFT matches a naive DFT") { - for (int n : {2, 4, 8, 16, 64, 256}) { - GIVEN("random complex input of length " + std::to_string(n)) { - std::vector re = noise(n, 100u + n), im = noise(n, 900u + n); - std::vector ref_re, ref_im; - naive_dft(re, im, ref_re, ref_im); - - std::vector fre = re, fim = im; - tap::tools::fft::transform(fre, fim, false); - - THEN("the outputs agree to within numerical tolerance") { - REQUIRE(max_abs_diff(fre, ref_re) < 1e-9); - REQUIRE(max_abs_diff(fim, ref_im) < 1e-9); - } - } - } -} - -SCENARIO("forward followed by inverse reconstructs the input") { - const int n = 128; - std::vector re = noise(n, 7u), im = noise(n, 8u); - const std::vector re0 = re, im0 = im; - - tap::tools::fft::transform(re, im, false); - tap::tools::fft::transform(re, im, true); - - THEN("the round trip returns the original (inverse divides by N)") { - REQUIRE(max_abs_diff(re, re0) < 1e-9); - REQUIRE(max_abs_diff(im, im0) < 1e-9); - } -} - -SCENARIO("a unit impulse transforms to a flat unit spectrum") { - const int n = 16; - std::vector re(n, 0.0), im(n, 0.0); - re[0] = 1.0; - tap::tools::fft::transform(re, im, false); - bool flat = true; - for (int k = 0; k < n; ++k) { - if (std::abs(re[k] - 1.0) > 1e-12 || std::abs(im[k]) > 1e-12) { - flat = false; - } - } - REQUIRE(flat); -} - -SCENARIO("a real cosine lands on the two symmetric bins") { - const int n = 32; - const int freq = 3; - std::vector re(n), im(n, 0.0); - for (int t = 0; t < n; ++t) { - re[t] = std::cos(2.0 * k_pi * freq * t / n); - } - tap::tools::fft::transform(re, im, false); - - THEN("energy is N/2 at bin freq and N-freq, ~0 elsewhere") { - REQUIRE(std::abs(re[freq] - n / 2.0) < 1e-9); - REQUIRE(std::abs(re[n - freq] - n / 2.0) < 1e-9); - for (int k = 0; k < n; ++k) { - if (k != freq && k != n - freq) { - REQUIRE(std::abs(re[k]) < 1e-9); - } - } - } -} diff --git a/tests/spectra_test.cpp b/tests/spectra_test.cpp index 3da1e71..e5484ce 100644 --- a/tests/spectra_test.cpp +++ b/tests/spectra_test.cpp @@ -9,7 +9,7 @@ #include #include -#include +#include #include namespace { @@ -34,14 +34,16 @@ namespace { return v; } - // Peak magnitude bin of an N-sample slice of `sig` starting at `off` (rectangular window). + // Peak magnitude bin of an N-sample slice of `sig` starting at `off` (rectangular window), + // via the shared DspTap real FFT's packed spectrum (DC in slot 0, Nyquist in slot 1). int peak_bin(const std::vector& sig, int off, int N) { - std::vector re(sig.begin() + off, sig.begin() + off + N), im(N, 0.0); - tap::tools::fft::transform(re, im, false); + std::vector a(sig.begin() + off, sig.begin() + off + N); + tap::dsp::real_fft fft(static_cast(N)); + fft.forward_inplace(a.data()); int best = 0; double bestmag = -1.0; for (int k = 1; k <= N / 2; ++k) { - const double mag = re[k] * re[k] + im[k] * im[k]; + const double mag = (k == N / 2) ? (a[1] * a[1]) : (a[2 * k] * a[2 * k] + a[2 * k + 1] * a[2 * k + 1]); if (mag > bestmag) { bestmag = mag; best = k; diff --git a/tools/capi/CMakeLists.txt b/tools/capi/CMakeLists.txt index 7ee2e42..aaa4b72 100644 --- a/tools/capi/CMakeLists.txt +++ b/tools/capi/CMakeLists.txt @@ -14,6 +14,13 @@ project(taptools_capi CXX) add_library(taptools_capi SHARED taptools_capi.cpp) target_include_directories(taptools_capi PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../include) +# conv_engine now uses the shared DspTap real FFT (tap::dsp). When this builds as part of the +# kernel project the target already exists; standalone, pull the submodule in directly. +if (NOT TARGET tap::dsp) + add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/../../submodules/dsptap ${CMAKE_CURRENT_BINARY_DIR}/dsptap) +endif () +target_link_libraries(taptools_capi PRIVATE tap::dsp) + set_target_properties(taptools_capi PROPERTIES CXX_STANDARD 20 CXX_STANDARD_REQUIRED ON