diff --git a/CMakeLists.txt b/CMakeLists.txt index a0569ef..a4ecebb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.15) +cmake_minimum_required(VERSION 3.5.0) if (NOT DEFINED GUI_BASE_DIR) if (DEFINED ENV{GUI_BASE_DIR}) set(GUI_BASE_DIR $ENV{GUI_BASE_DIR}) diff --git a/DependencyTests/CMakeLists.txt b/DependencyTests/CMakeLists.txt index 167c41b..7df81c3 100644 --- a/DependencyTests/CMakeLists.txt +++ b/DependencyTests/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.15) +cmake_minimum_required(VERSION 3.22) project(OpenEphysFFTWDependencySmokeTest LANGUAGES C) if(NOT DEFINED FFTW_ROOT) diff --git a/README.md b/README.md index f387150..015c3ae 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,12 @@ This is a small library to wrap FFTW functionality useful for Open Ephys. It is currently used by the [Phase Calculator](https://github.com/tne-lab/phase-calculator) and Real-Time Coherence plugins (under development). +`OpenEphysFFTWBatch.h` provides a lightweight, JUCE-independent batch API for +aligned single- and double-precision real-to-complex and complex-to-real +transforms. Batches are transform-major, with contiguous rows. Inverse +transforms follow FFTW's convention and are not normalized. Create plans off +real-time threads and reuse them. + ## Dependency OpenEphysFFTW uses complete non-threaded, shared double- and single-precision @@ -24,8 +30,8 @@ previously staged platform directory: cmake -S . -B Build -DFFTW_ROOT=/path/to/libs/linux ``` -Automatic staging requires CMake 3.24 or newer. Builds using `FFTW_ROOT` retain -the project's CMake 3.15 minimum. +Automatic staging requires CMake 3.19 or newer. Builds using `FFTW_ROOT` retain +the project's existing CMake 3.5 minimum. The bundles intentionally omit the separate long-double, MPI, OpenMP, and threaded FFTW libraries. The `fftw3` and `fftw3f` libraries themselves are not diff --git a/Source/OpenEphysFFTW.cpp b/Source/OpenEphysFFTW.cpp index 4a1fb18..fb58c45 100644 --- a/Source/OpenEphysFFTW.cpp +++ b/Source/OpenEphysFFTW.cpp @@ -353,4 +353,4 @@ void FFTWTransformableArray::hilbert() fftReal(); freqDomainHilbert(); ifft(); -} \ No newline at end of file +} diff --git a/Source/OpenEphysFFTW.h b/Source/OpenEphysFFTW.h index 36d913d..8331fab 100644 --- a/Source/OpenEphysFFTW.h +++ b/Source/OpenEphysFFTW.h @@ -32,6 +32,8 @@ transform library #include #include +#include "OpenEphysFFTWBatch.h" + // forward-declare: struct fftw_plan_s; @@ -174,4 +176,4 @@ class FFTWTransformableArrayUsing : public FFTWTransformableArray {} }; -#endif // OEP_FFTW_H_INCLUDED \ No newline at end of file +#endif // OEP_FFTW_H_INCLUDED diff --git a/Source/OpenEphysFFTWBatch.cpp b/Source/OpenEphysFFTWBatch.cpp new file mode 100644 index 0000000..046f9b4 --- /dev/null +++ b/Source/OpenEphysFFTWBatch.cpp @@ -0,0 +1,422 @@ +/* +------------------------------------------------------------------ + +This file is part of a library for the Open Ephys GUI +Copyright (C) 2026 Open Ephys + +------------------------------------------------------------------ + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . + +*/ + +#include "OpenEphysFFTWBatch.h" + +#include + +#include +#include +#include + +namespace +{ +std::mutex plannerMutex; + +std::size_t checkedElementCount (int elementsPerTransform, int transformCount) +{ + if (elementsPerTransform <= 0 || transformCount <= 0) + throw std::invalid_argument ("FFTW batch dimensions must be positive"); + + const auto elements = static_cast (elementsPerTransform); + const auto transforms = static_cast (transformCount); + if (transforms > std::numeric_limits::max() / elements) + throw std::length_error ("FFTW batch allocation is too large"); + return elements * transforms; +} + +bool validTransformIndex (int transformIndex, int transformCount) +{ + return transformIndex >= 0 && transformIndex < transformCount; +} +} // namespace + +struct FFTWRealToComplexBatchFloat::Impl +{ + Impl (int length, int count, unsigned int flags) + : transformLength (length), transformCount (count), binCount (length > 0 ? length / 2 + 1 : 0), input (fftwf_alloc_real (checkedElementCount (length, count))), output (fftwf_alloc_complex (checkedElementCount (binCount, count))) + { + if (input == nullptr || output == nullptr) + { + fftwf_free (output); + fftwf_free (input); + throw std::bad_alloc(); + } + + const int dimensions[] { transformLength }; + { + const std::lock_guard lock (plannerMutex); + plan = fftwf_plan_many_dft_r2c (1, dimensions, transformCount, input, nullptr, 1, transformLength, output, nullptr, 1, binCount, flags); + } + if (plan == nullptr) + { + fftwf_free (output); + fftwf_free (input); + throw std::runtime_error ( + "Unable to create FFTW single-precision batch plan"); + } + } + + ~Impl() + { + if (plan != nullptr) + { + const std::lock_guard lock (plannerMutex); + fftwf_destroy_plan (plan); + } + fftwf_free (output); + fftwf_free (input); + } + + const int transformLength; + const int transformCount; + const int binCount; + float* input = nullptr; + fftwf_complex* output = nullptr; + fftwf_plan plan = nullptr; +}; + +FFTWRealToComplexBatchFloat::FFTWRealToComplexBatchFloat (int transformLength, + int transformCount, + unsigned int flags) + : impl (new Impl (transformLength, transformCount, flags)) {} + +FFTWRealToComplexBatchFloat::~FFTWRealToComplexBatchFloat() = default; + +float* FFTWRealToComplexBatchFloat::getInputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->input + transformIndex * impl->transformLength + : nullptr; +} + +const float* + FFTWRealToComplexBatchFloat::getInputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->input + transformIndex * impl->transformLength + : nullptr; +} + +std::complex* + FFTWRealToComplexBatchFloat::getOutputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->output) + transformIndex * impl->binCount + : nullptr; +} + +const std::complex* + FFTWRealToComplexBatchFloat::getOutputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->output) + transformIndex * impl->binCount + : nullptr; +} + +int FFTWRealToComplexBatchFloat::getTransformLength() const noexcept +{ + return impl->transformLength; +} +int FFTWRealToComplexBatchFloat::getTransformCount() const noexcept +{ + return impl->transformCount; +} +int FFTWRealToComplexBatchFloat::getBinCount() const noexcept +{ + return impl->binCount; +} + +void FFTWRealToComplexBatchFloat::execute() { fftwf_execute (impl->plan); } + +struct FFTWRealToComplexBatchDouble::Impl +{ + Impl (int length, int count, unsigned int flags) + : transformLength (length), transformCount (count), binCount (length > 0 ? length / 2 + 1 : 0), input (fftw_alloc_real (checkedElementCount (length, count))), output (fftw_alloc_complex (checkedElementCount (binCount, count))) + { + if (input == nullptr || output == nullptr) + { + fftw_free (output); + fftw_free (input); + throw std::bad_alloc(); + } + + const int dimensions[] { transformLength }; + { + const std::lock_guard lock (plannerMutex); + plan = fftw_plan_many_dft_r2c (1, dimensions, transformCount, input, nullptr, 1, transformLength, output, nullptr, 1, binCount, flags); + } + if (plan == nullptr) + { + fftw_free (output); + fftw_free (input); + throw std::runtime_error ( + "Unable to create FFTW double-precision batch plan"); + } + } + + ~Impl() + { + if (plan != nullptr) + { + const std::lock_guard lock (plannerMutex); + fftw_destroy_plan (plan); + } + fftw_free (output); + fftw_free (input); + } + + const int transformLength; + const int transformCount; + const int binCount; + double* input = nullptr; + fftw_complex* output = nullptr; + fftw_plan plan = nullptr; +}; + +FFTWRealToComplexBatchDouble::FFTWRealToComplexBatchDouble (int transformLength, + int transformCount, + unsigned int flags) + : impl (new Impl (transformLength, transformCount, flags)) {} + +FFTWRealToComplexBatchDouble::~FFTWRealToComplexBatchDouble() = default; + +double* FFTWRealToComplexBatchDouble::getInputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->input + transformIndex * impl->transformLength + : nullptr; +} + +const double* + FFTWRealToComplexBatchDouble::getInputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->input + transformIndex * impl->transformLength + : nullptr; +} + +std::complex* + FFTWRealToComplexBatchDouble::getOutputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->output) + transformIndex * impl->binCount + : nullptr; +} + +const std::complex* + FFTWRealToComplexBatchDouble::getOutputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->output) + transformIndex * impl->binCount + : nullptr; +} + +int FFTWRealToComplexBatchDouble::getTransformLength() const noexcept +{ + return impl->transformLength; +} +int FFTWRealToComplexBatchDouble::getTransformCount() const noexcept +{ + return impl->transformCount; +} +int FFTWRealToComplexBatchDouble::getBinCount() const noexcept +{ + return impl->binCount; +} + +void FFTWRealToComplexBatchDouble::execute() { fftw_execute (impl->plan); } + +struct FFTWComplexToRealBatchFloat::Impl +{ + Impl (int length, int count, unsigned int flags) + : transformLength (length), transformCount (count), binCount (length > 0 ? length / 2 + 1 : 0), input (fftwf_alloc_complex (checkedElementCount (binCount, count))), output (fftwf_alloc_real (checkedElementCount (length, count))) + { + if (input == nullptr || output == nullptr) + { + fftwf_free (output); + fftwf_free (input); + throw std::bad_alloc(); + } + + const int dimensions[] { transformLength }; + { + const std::lock_guard lock (plannerMutex); + plan = fftwf_plan_many_dft_c2r (1, dimensions, transformCount, input, nullptr, 1, binCount, output, nullptr, 1, transformLength, flags); + } + if (plan == nullptr) + { + fftwf_free (output); + fftwf_free (input); + throw std::runtime_error ( + "Unable to create FFTW single-precision inverse batch plan"); + } + } + + ~Impl() + { + if (plan != nullptr) + { + const std::lock_guard lock (plannerMutex); + fftwf_destroy_plan (plan); + } + fftwf_free (output); + fftwf_free (input); + } + + const int transformLength; + const int transformCount; + const int binCount; + fftwf_complex* input = nullptr; + float* output = nullptr; + fftwf_plan plan = nullptr; +}; + +FFTWComplexToRealBatchFloat::FFTWComplexToRealBatchFloat (int transformLength, + int transformCount, + unsigned int flags) + : impl (new Impl (transformLength, transformCount, flags)) {} + +FFTWComplexToRealBatchFloat::~FFTWComplexToRealBatchFloat() = default; + +std::complex* FFTWComplexToRealBatchFloat::getInputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->input) + transformIndex * impl->binCount + : nullptr; +} + +const std::complex* + FFTWComplexToRealBatchFloat::getInputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->input) + transformIndex * impl->binCount + : nullptr; +} + +float* FFTWComplexToRealBatchFloat::getOutputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->output + transformIndex * impl->transformLength + : nullptr; +} + +const float* FFTWComplexToRealBatchFloat::getOutputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->output + transformIndex * impl->transformLength + : nullptr; +} + +int FFTWComplexToRealBatchFloat::getTransformLength() const noexcept { return impl->transformLength; } +int FFTWComplexToRealBatchFloat::getTransformCount() const noexcept { return impl->transformCount; } +int FFTWComplexToRealBatchFloat::getBinCount() const noexcept { return impl->binCount; } + +void FFTWComplexToRealBatchFloat::execute() { fftwf_execute (impl->plan); } + +struct FFTWComplexToRealBatchDouble::Impl +{ + Impl (int length, int count, unsigned int flags) + : transformLength (length), transformCount (count), binCount (length > 0 ? length / 2 + 1 : 0), input (fftw_alloc_complex (checkedElementCount (binCount, count))), output (fftw_alloc_real (checkedElementCount (length, count))) + { + if (input == nullptr || output == nullptr) + { + fftw_free (output); + fftw_free (input); + throw std::bad_alloc(); + } + + const int dimensions[] { transformLength }; + { + const std::lock_guard lock (plannerMutex); + plan = fftw_plan_many_dft_c2r (1, dimensions, transformCount, input, nullptr, 1, binCount, output, nullptr, 1, transformLength, flags); + } + if (plan == nullptr) + { + fftw_free (output); + fftw_free (input); + throw std::runtime_error ( + "Unable to create FFTW double-precision inverse batch plan"); + } + } + + ~Impl() + { + if (plan != nullptr) + { + const std::lock_guard lock (plannerMutex); + fftw_destroy_plan (plan); + } + fftw_free (output); + fftw_free (input); + } + + const int transformLength; + const int transformCount; + const int binCount; + fftw_complex* input = nullptr; + double* output = nullptr; + fftw_plan plan = nullptr; +}; + +FFTWComplexToRealBatchDouble::FFTWComplexToRealBatchDouble (int transformLength, + int transformCount, + unsigned int flags) + : impl (new Impl (transformLength, transformCount, flags)) {} + +FFTWComplexToRealBatchDouble::~FFTWComplexToRealBatchDouble() = default; + +std::complex* FFTWComplexToRealBatchDouble::getInputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->input) + transformIndex * impl->binCount + : nullptr; +} + +const std::complex* + FFTWComplexToRealBatchDouble::getInputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? reinterpret_cast*> (impl->input) + transformIndex * impl->binCount + : nullptr; +} + +double* FFTWComplexToRealBatchDouble::getOutputPointer (int transformIndex) +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->output + transformIndex * impl->transformLength + : nullptr; +} + +const double* FFTWComplexToRealBatchDouble::getOutputPointer (int transformIndex) const +{ + return validTransformIndex (transformIndex, impl->transformCount) + ? impl->output + transformIndex * impl->transformLength + : nullptr; +} + +int FFTWComplexToRealBatchDouble::getTransformLength() const noexcept { return impl->transformLength; } +int FFTWComplexToRealBatchDouble::getTransformCount() const noexcept { return impl->transformCount; } +int FFTWComplexToRealBatchDouble::getBinCount() const noexcept { return impl->binCount; } + +void FFTWComplexToRealBatchDouble::execute() { fftw_execute (impl->plan); } diff --git a/Source/OpenEphysFFTWBatch.h b/Source/OpenEphysFFTWBatch.h new file mode 100644 index 0000000..f404f48 --- /dev/null +++ b/Source/OpenEphysFFTWBatch.h @@ -0,0 +1,159 @@ +/* +------------------------------------------------------------------ + +This file is part of a library for the Open Ephys GUI +Copyright (C) 2026 Open Ephys + +------------------------------------------------------------------ + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . + +*/ + +#ifndef OEP_FFTW_BATCH_H_INCLUDED +#define OEP_FFTW_BATCH_H_INCLUDED + +#include +#include +#include + +/** + Owns aligned storage and a reusable single-precision real-to-complex batch. + + Each transform occupies one contiguous row of transformLength real inputs + and produces transformLength / 2 + 1 contiguous complex outputs. Planning + and allocation occur in the constructor; execute() performs no allocation. + The internal FFTW planning strategy is not part of this API contract. + Planning is serialized across batch instances. Concurrent execution is safe + for distinct instances, but callers must not execute one instance twice at + the same time because its input and output storage are shared. +*/ +class COMMON_LIB FFTWRealToComplexBatchFloat +{ +public: + FFTWRealToComplexBatchFloat (int transformLength, + int transformCount, + unsigned int flags = 0U /* FFTW_MEASURE */); + ~FFTWRealToComplexBatchFloat(); + + float* getInputPointer (int transformIndex = 0); + const float* getInputPointer (int transformIndex = 0) const; + std::complex* getOutputPointer (int transformIndex = 0); + const std::complex* getOutputPointer (int transformIndex = 0) const; + + int getTransformLength() const noexcept; + int getTransformCount() const noexcept; + int getBinCount() const noexcept; + + void execute(); + +private: + struct Impl; + std::unique_ptr impl; + + FFTWRealToComplexBatchFloat (const FFTWRealToComplexBatchFloat&) = delete; + FFTWRealToComplexBatchFloat& operator= (const FFTWRealToComplexBatchFloat&) = delete; +}; + +/** Double-precision counterpart to FFTWRealToComplexBatchFloat. */ +class COMMON_LIB FFTWRealToComplexBatchDouble +{ +public: + FFTWRealToComplexBatchDouble (int transformLength, + int transformCount, + unsigned int flags = 0U /* FFTW_MEASURE */); + ~FFTWRealToComplexBatchDouble(); + + double* getInputPointer (int transformIndex = 0); + const double* getInputPointer (int transformIndex = 0) const; + std::complex* getOutputPointer (int transformIndex = 0); + const std::complex* getOutputPointer (int transformIndex = 0) const; + + int getTransformLength() const noexcept; + int getTransformCount() const noexcept; + int getBinCount() const noexcept; + + void execute(); + +private: + struct Impl; + std::unique_ptr impl; + + FFTWRealToComplexBatchDouble (const FFTWRealToComplexBatchDouble&) = delete; + FFTWRealToComplexBatchDouble& operator= (const FFTWRealToComplexBatchDouble&) = delete; +}; + +/** + Reusable single-precision complex-to-real batch. + + Input rows contain transformLength / 2 + 1 complex values in FFTW's + Hermitian half-spectrum layout; output rows contain transformLength real + values. Like FFTW itself, execute() does not normalize the inverse transform. +*/ +class COMMON_LIB FFTWComplexToRealBatchFloat +{ +public: + FFTWComplexToRealBatchFloat (int transformLength, + int transformCount, + unsigned int flags = 0U /* FFTW_MEASURE */); + ~FFTWComplexToRealBatchFloat(); + + std::complex* getInputPointer (int transformIndex = 0); + const std::complex* getInputPointer (int transformIndex = 0) const; + float* getOutputPointer (int transformIndex = 0); + const float* getOutputPointer (int transformIndex = 0) const; + + int getTransformLength() const noexcept; + int getTransformCount() const noexcept; + int getBinCount() const noexcept; + + void execute(); + +private: + struct Impl; + std::unique_ptr impl; + + FFTWComplexToRealBatchFloat (const FFTWComplexToRealBatchFloat&) = delete; + FFTWComplexToRealBatchFloat& operator= (const FFTWComplexToRealBatchFloat&) = delete; +}; + +/** Double-precision counterpart to FFTWComplexToRealBatchFloat. */ +class COMMON_LIB FFTWComplexToRealBatchDouble +{ +public: + FFTWComplexToRealBatchDouble (int transformLength, + int transformCount, + unsigned int flags = 0U /* FFTW_MEASURE */); + ~FFTWComplexToRealBatchDouble(); + + std::complex* getInputPointer (int transformIndex = 0); + const std::complex* getInputPointer (int transformIndex = 0) const; + double* getOutputPointer (int transformIndex = 0); + const double* getOutputPointer (int transformIndex = 0) const; + + int getTransformLength() const noexcept; + int getTransformCount() const noexcept; + int getBinCount() const noexcept; + + void execute(); + +private: + struct Impl; + std::unique_ptr impl; + + FFTWComplexToRealBatchDouble (const FFTWComplexToRealBatchDouble&) = delete; + FFTWComplexToRealBatchDouble& operator= (const FFTWComplexToRealBatchDouble&) = delete; +}; + +#endif // OEP_FFTW_BATCH_H_INCLUDED diff --git a/tools/PrepareFFTWBundle.cmake b/tools/PrepareFFTWBundle.cmake index 918813d..19c7c43 100644 --- a/tools/PrepareFFTWBundle.cmake +++ b/tools/PrepareFFTWBundle.cmake @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.24) +cmake_minimum_required(VERSION 3.19) if(NOT DEFINED PLATFORM OR NOT PLATFORM MATCHES "^(linux|macos|windows)$") message(FATAL_ERROR "Set PLATFORM to linux, macos, or windows")