From 2d1fad7b33ea1ad6f392d48815efb1bead53cf2a Mon Sep 17 00:00:00 2001 From: slebedev Date: Fri, 22 Sep 2023 17:02:47 +0200 Subject: [PATCH] Itarative radix-2 FFT algorithm implementation * using but-reversal permutation. * FFT block contains 3 template parameters: data type and algorithm type. * Add test of different FFT algorithms. --- CMakeLists.txt | 3 +- bench/bm_fft.cpp | 129 +++---- test/blocklib/algorithm/fft/fft.hpp | 112 ++++++ test/blocklib/algorithm/fft/fft_common.hpp | 34 ++ test/blocklib/algorithm/fft/fft_types.hpp | 37 ++ test/blocklib/algorithm/fft/fftw.hpp | 215 +++++++++++ .../{core => algorithm}/fft/window.hpp | 23 +- test/blocklib/core/fft/fft.hpp | 353 ++++++------------ test/qa_fft.cpp | 257 ++++++++----- 9 files changed, 734 insertions(+), 429 deletions(-) create mode 100644 test/blocklib/algorithm/fft/fft.hpp create mode 100644 test/blocklib/algorithm/fft/fft_common.hpp create mode 100644 test/blocklib/algorithm/fft/fft_types.hpp create mode 100644 test/blocklib/algorithm/fft/fftw.hpp rename test/blocklib/{core => algorithm}/fft/window.hpp (91%) diff --git a/CMakeLists.txt b/CMakeLists.txt index e8f84a8f9..75c999ad5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -85,11 +85,12 @@ target_include_directories(pmtv INTERFACE ${pmt_SOURCE_DIR}/include/) # FFTW3 is build 2 times for float and double precisions SET(FFTW_PREFIX ${CMAKE_CURRENT_BINARY_DIR}/fftw) -SET(FFTW_CFLAGS "${CFLAGS} -w") if (EMSCRIPTEN) + SET(FFTW_CFLAGS "${CFLAGS} -w") SET(FFTW_CONFIG cd ${FFTW_PREFIX}/src/ && emconfigure ./configure --enable-silent-rules --quiet --disable-fortran --prefix=${FFTW_PREFIX}/install) SET(FFTW_BUILD emmake make -j CFLAGS=${FFTW_CFLAGS} --silent V=0 && emmake make install --silent V=0 && emmake make clean --silent V=0) else () + SET(FFTW_CFLAGS "${CFLAGS} -w -O3 -march=native -mtune=native") SET(FFTW_CONFIG ${FFTW_PREFIX}/src/configure --enable-silent-rules --quiet --disable-fortran --prefix=${FFTW_PREFIX}/install) SET(FFTW_BUILD make -j CFLAGS=${FFTW_CFLAGS} --silent V=0 && make install --silent V=0 && make clean --silent V=0) endif () diff --git a/bench/bm_fft.cpp b/bench/bm_fft.cpp index 5cb107d31..2cab1115d 100644 --- a/bench/bm_fft.cpp +++ b/bench/bm_fft.cpp @@ -1,53 +1,19 @@ -#include "benchmark.hpp" - +#include "../test/blocklib/algorithm/fft/fft.hpp" +#include "../test/blocklib/algorithm/fft/fftw.hpp" #include "../test/blocklib/core/fft/fft.hpp" - +#include "benchmark.hpp" #include #include -/// This custom implementation of FFT ("compute_fft_v1" and "computeFftV1") is done only for performance comparison with default FFTW implementation. - -/** - * Real-valued fast fourier transform algorithms - * H.V. Sorensen, D.L. Jones, M.T. Heideman, C.S. Burrus (1987), - * in: IEEE Trans on Acoustics, Speech, & Signal Processing, 35 - */ - -template - requires gr::blocks::fft::ComplexType -void -computeFftV1(std::vector &signal) { - const std::size_t N{ signal.size() }; - for (std::size_t j = 0, rev = 0; j < N; j++) { - if (j < rev) std::swap(signal[j], signal[rev]); - auto maskLen = static_cast(std::countr_zero(j + 1) + 1); - rev ^= N - (N >> maskLen); - } - - for (std::size_t s = 2; s <= N; s *= 2) { - const std::size_t m{ s / 2 }; - const T w{ exp(T(0., static_cast(-2. * std::numbers::pi) / static_cast(s))) }; - for (std::size_t k = 0; k < N; k += s) { - T wk{ 1., 0. }; - for (std::size_t j = 0; j < m; j++) { - const T t{ wk * signal[k + j + m] }; - const T u{ signal[k + j] }; - signal[k + j] = u + t; - signal[k + j + m] = u - t; - wk *= w; - } - } - } -} - +/// This custom implementation of FFT is done only for performance comparison with default FFTW implementation. /** * Fast Fourier-Transform according to Cooley–Tukey * reference: https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#Pseudocode */ template - requires gr::blocks::fft::ComplexType + requires gr::algorithm::ComplexType void -computeFftV2(std::vector &signal) { +computeFFFTCooleyTukey(std::vector &signal) { const std::size_t N{ signal.size() }; if (N == 1) return; @@ -59,10 +25,10 @@ computeFftV2(std::vector &signal) { odd[i] = signal[2 * i + 1]; } - computeFftV2(even); - computeFftV2(odd); + computeFFFTCooleyTukey(even); + computeFFFTCooleyTukey(odd); - const typename T::value_type wn{ static_cast(2. * std::numbers::pi) / static_cast(N) }; + const typename T::value_type wn{ static_cast(2. * std::numbers::pi_v) / static_cast(N) }; for (std::size_t i = 0; i < N / 2; i++) { const T wkn(std::cos(wn * static_cast(i)), std::sin(wn * static_cast(i))); signal[i] = even[i] + wkn * odd[i]; @@ -70,24 +36,12 @@ computeFftV2(std::vector &signal) { } } -template - requires gr::blocks::fft::ComplexType -std::vector -computeMagnitudeSpectrum(std::vector &fftSignal) { - const std::size_t N{ fftSignal.size() }; - std::vector magnitudeSpectrum(N / 2); - for (std::size_t i = 0; i < N / 2; i++) { - magnitudeSpectrum[i] = std::hypot(fftSignal[i].real(), fftSignal[i].imag()) * static_cast(2.) / static_cast(N); - } - return magnitudeSpectrum; -} - template std::vector generateSinSample(std::size_t N, double sampleRate, double frequency, double amplitude) { std::vector signal(N); for (std::size_t i = 0; i < N; i++) { - if constexpr (gr::blocks::fft::ComplexType) { + if constexpr (gr::algorithm::ComplexType) { signal[i] = { static_cast(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast(i) / sampleRate)), 0. }; } else { signal[i] = static_cast(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast(i) / sampleRate)); @@ -98,7 +52,7 @@ generateSinSample(std::size_t N, double sampleRate, double frequency, double amp template void -testFft(bool withMagSpectrum) { +testFFT() { using namespace benchmark; using namespace boost::ut::reflection; @@ -110,50 +64,67 @@ testFft(bool withMagSpectrum) { static_assert(std::has_single_bit(N)); - std::vector signal = generateSinSample(N, sampleRate, frequency, amplitude); - std::string nameOpt = withMagSpectrum ? "fft+mag" : "fft"; + std::vector signal = generateSinSample(N, sampleRate, frequency, amplitude); { - gr::blocks::fft::fft fft1{}; + gr::blocks::fft::FFT> fft1{}; std::ignore = fft1.settings().set({ { "fftSize", N } }); std::ignore = fft1.settings().apply_staged_parameters(); fft1.inputHistory.push_back_bulk(signal.begin(), signal.end()); - ::benchmark::benchmark(fmt::format("{} - {} fftw", nameOpt, type_name())) = [&fft1, &withMagSpectrum] { + ::benchmark::benchmark(fmt::format("{} - fftw", type_name())) = [&fft1] { + fft1.prepareInput(); + fft1.computeFFT(); + }; + } + { + gr::blocks::fft::FFT> fft1; + std::ignore = fft1.settings().set({ { "fftSize", N } }); + std::ignore = fft1.settings().apply_staged_parameters(); + fft1.inputHistory.push_back_bulk(signal.begin(), signal.end()); + ::benchmark::benchmark(fmt::format("{} - fft", type_name())) = [&fft1] { fft1.prepareInput(); - fft1.computeFft(); - if (withMagSpectrum) fft1.computeMagnitudeSpectrum(); + fft1.computeFFT(); }; } if constexpr (gr::blocks::fft::ComplexType) { - ::benchmark::benchmark(fmt::format("{} - {} fft_v1", nameOpt, type_name())) = [&signal, &withMagSpectrum] { + ::benchmark::benchmark(fmt::format("{} - fftCT", type_name())) = [&signal] { auto signalCopy = signal; - computeFftV1(signalCopy); - if (withMagSpectrum) [[maybe_unused]] - auto magnitudeSpectrum = computeMagnitudeSpectrum(signalCopy); + computeFFFTCooleyTukey(signalCopy); }; + } - ::benchmark::benchmark(fmt::format("{} - {} fft_v2", nameOpt, type_name())) = [&signal, &withMagSpectrum] { - auto signalCopy = signal; - computeFftV2(signalCopy); - if (withMagSpectrum) [[maybe_unused]] - auto magnitudeSpectrum = computeMagnitudeSpectrum(signalCopy); + { + gr::blocks::fft::FFT> fft1{}; + std::ignore = fft1.settings().set({ { "fftSize", N } }); + std::ignore = fft1.settings().apply_staged_parameters(); + fft1.inputHistory.push_back_bulk(signal.begin(), signal.end()); + fft1.prepareInput(); + fft1.computeFFT(); + + ::benchmark::benchmark(fmt::format("{} - magnitude", type_name())) = [&fft1] { fft1.computeMagnitudeSpectrum(); }; + ::benchmark::benchmark(fmt::format("{} - phase", type_name())) = [&fft1] { fft1.computePhaseSpectrum(); }; + + double sum{ 0. }; + fft1.computeMagnitudeSpectrum(); + fft1.computePhaseSpectrum(); + ::benchmark::benchmark(fmt::format("{} - dataset", type_name())) = [&fft1, &sum] { + const auto ds1 = fft1.createDataset(); + sum += static_cast(ds1.signal_values[0]); }; } + + ::benchmark::results::add_separator(); } inline const boost::ut::suite _fft_bm_tests = [] { std::tuple, std::complex> complexTypesToTest{}; std::tuple realTypesToTest{}; - std::apply([](TArgs... /*args*/) { (testFft(false), ...); }, complexTypesToTest); - benchmark::results::add_separator(); - std::apply([](TArgs... /*args*/) { (testFft(true), ...); }, complexTypesToTest); - benchmark::results::add_separator(); - std::apply([](TArgs... /*args*/) { (testFft(false), ...); }, realTypesToTest); - benchmark::results::add_separator(); - std::apply([](TArgs... /*args*/) { (testFft(true), ...); }, realTypesToTest); + std::apply([](TArgs... /*args*/) { (testFFT(), ...); }, complexTypesToTest); + // benchmark::results::add_separator(); + std::apply([](TArgs... /*args*/) { (testFFT(), ...); }, realTypesToTest); }; int diff --git a/test/blocklib/algorithm/fft/fft.hpp b/test/blocklib/algorithm/fft/fft.hpp new file mode 100644 index 000000000..cb9d5485b --- /dev/null +++ b/test/blocklib/algorithm/fft/fft.hpp @@ -0,0 +1,112 @@ +#ifndef GRAPH_PROTOTYPE_ALGORITHM_FFT_HPP +#define GRAPH_PROTOTYPE_ALGORITHM_FFT_HPP + +#include "fft_types.hpp" +#include "node.hpp" +#include "window.hpp" +#include + +namespace gr::algorithm { + +template + requires(ComplexType || std::floating_point) +struct FFT { + using PrecisionType = typename FFTInternal::PrecisionType; + using OutDataType = typename FFTInternal::OutDataType; + + std::vector twiddleFactors{}; + std::size_t fftSize{ 0 }; + + FFT() = default; + FFT(const FFT &rhs) = delete; + FFT(FFT &&rhs) noexcept = delete; + FFT & + operator=(const FFT &rhs) + = delete; + FFT & + operator=(FFT &&rhs) noexcept + = delete; + + ~FFT() = default; + + void + initAll() { + precomputeTwiddleFactors(); + } + + std::vector + computeFFT(const std::vector &in) noexcept { + std::vector out(in.size()); + computeFFT(in, out); + return out; + } + + void + computeFFT(const std::vector &in, std::vector &out) noexcept { + if (!std::has_single_bit(in.size())) { + throw std::runtime_error(fmt::format("Input data must have 2^N samples, input size: ", in.size())); + } + if (fftSize != in.size()) { + fftSize = in.size(); + initAll(); + } + + // For the moment no optimization for real data inputs, just create complex with zero imaginary value. + if constexpr (!ComplexType) { + std::transform(in.begin(), in.end(), out.begin(), [](const auto c) { return OutDataType(c, 0.); }); + } else { + std::copy(in.begin(), in.end(), out.begin()); + } + + /** + * Real-valued fast fourier transform algorithms + * H.V. Sorensen, D.L. Jones, M.T. Heideman, C.S. Burrus (1987), + * in: IEEE Trans on Acoustics, Speech, & Signal Processing, 35 + */ + bitReversalPermutation(out); + + std::size_t omega_kCounter = 0; + for (std::size_t s = 2; s <= fftSize; s *= 2) { + const auto half_s = s / 2; + for (std::size_t k = 0; k < fftSize; k += s) { + for (std::size_t j = 0; j < half_s; j++) { + const auto t{ twiddleFactors[omega_kCounter++] * out[k + j + half_s] }; + const auto u{ out[k + j] }; + out[k + j] = u + t; + out[k + j + half_s] = u - t; + } + } + } + } + +private: + void + bitReversalPermutation(std::vector &vec) const noexcept { + for (std::size_t j = 0, rev = 0; j < fftSize; j++) { + if (j < rev) std::swap(vec[j], vec[rev]); + auto maskLen = static_cast(std::countr_zero(j + 1) + 1); + rev ^= fftSize - (fftSize >> maskLen); + } + } + + void + precomputeTwiddleFactors() { + twiddleFactors.clear(); + const auto minus2Pi = PrecisionType(-2. * std::numbers::pi); + for (std::size_t s = 2; s <= fftSize; s *= 2) { + const std::size_t m{ s / 2 }; + const OutDataType w{ std::exp(OutDataType(0., minus2Pi / static_cast(s))) }; + for (std::size_t k = 0; k < fftSize; k += s) { + OutDataType wk{ 1., 0. }; + for (std::size_t j = 0; j < m; j++) { + twiddleFactors.push_back(wk); + wk *= w; + } + } + } + } +}; + +} // namespace gr::algorithm + +#endif // GRAPH_PROTOTYPE_ALGORITHM_FFT_HPP diff --git a/test/blocklib/algorithm/fft/fft_common.hpp b/test/blocklib/algorithm/fft/fft_common.hpp new file mode 100644 index 000000000..c3a7cdc93 --- /dev/null +++ b/test/blocklib/algorithm/fft/fft_common.hpp @@ -0,0 +1,34 @@ +#ifndef GRAPH_PROTOTYPE_ALGORITHM_FFT_COMMON_HPP +#define GRAPH_PROTOTYPE_ALGORITHM_FFT_COMMON_HPP + +#include "fft_types.hpp" +#include +#include +#include + +namespace gr::algorithm { + +template +void +computeMagnitudeSpectrum(const std::vector &fftOut, std::vector &magnitudeSpectrum, std::size_t fftSize, bool outputInDb) { + if (fftOut.size() < magnitudeSpectrum.size()) { + throw std::invalid_argument(fmt::format("FFT vector size ({}) must be more or equal than magnitude vector size ({}).", fftOut.size(), magnitudeSpectrum.size())); + } + using PrecisionType = typename T::value_type; + std::transform(fftOut.begin(), std::next(fftOut.begin(), magnitudeSpectrum.size()), magnitudeSpectrum.begin(), [fftSize, outputInDb](const auto &c) { + const auto mag{ std::hypot(c.real(), c.imag()) * PrecisionType(2.0) / static_cast(fftSize) }; + return static_cast(outputInDb ? PrecisionType(20.) * std::log10(std::abs(mag)) : mag); + }); +} + +template +void +computePhaseSpectrum(const std::vector &fftOut, std::vector &phaseSpectrum) { + if (fftOut.size() < phaseSpectrum.size()) { + throw std::invalid_argument(fmt::format("FFT vector size ({}) must be more or equal than phase vector size ({}).", fftOut.size(), phaseSpectrum.size())); + } + std::transform(fftOut.begin(), std::next(fftOut.begin(), phaseSpectrum.size()), phaseSpectrum.begin(), [](const auto &c) { return static_cast(std::atan2(c.imag(), c.real())); }); +} + +} // namespace gr::algorithm +#endif // GRAPH_PROTOTYPE_ALGORITHM_FFT_COMMON_HPP diff --git a/test/blocklib/algorithm/fft/fft_types.hpp b/test/blocklib/algorithm/fft/fft_types.hpp new file mode 100644 index 000000000..2c134e6b0 --- /dev/null +++ b/test/blocklib/algorithm/fft/fft_types.hpp @@ -0,0 +1,37 @@ +#ifndef GRAPH_PROTOTYPE_ALGORITHM_FFT_TYPES_HPP +#define GRAPH_PROTOTYPE_ALGORITHM_FFT_TYPES_HPP + +#include "node.hpp" + +namespace gr::algorithm { +template +concept ComplexType = std::is_same_v> || std::is_same_v>; + +template +concept DoubleType = std::is_same_v> || std::is_same_v; + +// clang-format off +template struct MakeSigned { using type = signedT;}; +template struct MakeSigned { using type = std::make_signed_t; }; +template struct EvalOutputType { using type1 = evalU; using type2 = typename MakeSigned::type;}; +template struct EvalOutputType { using type1 = typename evalU::value_type; using type2 = evalU;}; + +template +struct FFTInternal { + using PrecisionType = float; + using InDataType = std::conditional_t, std::complex, float>; + using OutDataType = std::complex; +}; + +template +struct FFTInternal { + using PrecisionType = double; + using InDataType = std::conditional_t, std::complex, double>; + using OutDataType = std::complex; +}; + +template +using FFTOutputType = std::conditional_t, typename EvalOutputType::type1, typename EvalOutputType::type2>; + +} // namespace gr::algorithm +#endif // GRAPH_PROTOTYPE_ALGORITHM_FFT_TYPES_HPP diff --git a/test/blocklib/algorithm/fft/fftw.hpp b/test/blocklib/algorithm/fft/fftw.hpp new file mode 100644 index 000000000..5e6176bae --- /dev/null +++ b/test/blocklib/algorithm/fft/fftw.hpp @@ -0,0 +1,215 @@ +#ifndef GRAPH_PROTOTYPE_ALGORITHM_FFTW_HPP +#define GRAPH_PROTOTYPE_ALGORITHM_FFTW_HPP + +#include "dataset.hpp" +#include "fft_types.hpp" +#include "fftw3.h" +#include "history_buffer.hpp" +#include "node.hpp" +#include "window.hpp" + +namespace gr::algorithm { + +using namespace fair::graph; + +inline static std::mutex fftw_plan_mutex; + +template + requires(ComplexType || std::floating_point) +struct FFTw { +public: + // clang-format off + template + struct fftwImpl { + using PlanType = fftwf_plan; + using InAlgoDataType = std::conditional_t, fftwf_complex, float>; + using OutAlgoDataType = fftwf_complex; + using InUniquePtr = std::unique_ptr; + using OutUniquePtr = std::unique_ptr; + using PlanUniquePtr = std::unique_ptr, decltype([](PlanType ptr) { fftwf_destroy_plan(ptr); })>; + + static void execute(const PlanType p) { fftwf_execute(p); } + static void cleanup() { fftwf_cleanup(); } + static void * malloc(std::size_t n) { return fftwf_malloc(n);} + static PlanType plan(int p_n, InAlgoDataType *p_in, OutAlgoDataType *p_out, int p_sign, unsigned int p_flags) { + if constexpr (std::is_same_v) { + return fftwf_plan_dft_r2c_1d(p_n, p_in, p_out, p_flags); + } else { + return fftwf_plan_dft_1d(p_n, p_in, p_out, p_sign, p_flags); + } + } + static int importWisdomFromFilename(const std::string& path) {return fftwf_import_wisdom_from_filename(path.c_str());} + static int exportWisdomToFilename(const std::string& path) {return fftwf_export_wisdom_to_filename(path.c_str());} + static int importWisdomFromString(const std::string& str) {return fftwf_import_wisdom_from_string(str.c_str());} + static std::string exportWisdomToString() { + char* cstr = fftwf_export_wisdom_to_string(); + if (cstr == nullptr) return ""; + std::string str(cstr); + fftwf_free(cstr); + return str; + } + static void forgetWisdom() {fftwf_forget_wisdom();} + }; + + template + struct fftwImpl { + using PlanType = fftw_plan; + using InAlgoDataType = std::conditional_t, fftw_complex, double>; + using OutAlgoDataType = fftw_complex; + using InUniquePtr = std::unique_ptr; + using OutUniquePtr = std::unique_ptr; + using PlanUniquePtr = std::unique_ptr, decltype([](PlanType ptr) { fftw_destroy_plan(ptr); })>; + + static void execute(const PlanType p) { fftw_execute(p); } + static void cleanup() { fftw_cleanup(); } + static void * malloc(std::size_t n) { return fftw_malloc(n);} + static PlanType plan(int p_n, InAlgoDataType *p_in, OutAlgoDataType *p_out, int p_sign, unsigned int p_flags) { + if constexpr (std::is_same_v) { + return fftw_plan_dft_r2c_1d(p_n, p_in, p_out, p_flags); + } else { + return fftw_plan_dft_1d(p_n, p_in, p_out, p_sign, p_flags); + } + } + static int importWisdomFromFilename(const std::string& path) {return fftw_import_wisdom_from_filename(path.c_str());} + static int exportWisdomToFilename(const std::string& path) {return fftw_export_wisdom_to_filename(path.c_str());} + static int importWisdomFromString(const std::string& str) {return fftw_import_wisdom_from_string(str.c_str());} + static std::string exportWisdomToString() { + char* cstr = fftw_export_wisdom_to_string(); + if (cstr == nullptr) return ""; + std::string str(cstr); + fftw_free(cstr); + return str; + } + static void forgetWisdom() {fftw_forget_wisdom();} + }; + + // clang-format on + + using OutDataType = typename FFTInternal::OutDataType; + using InAlgoDataType = typename fftwImpl::InAlgoDataType; + using OutAlgoDataType = typename fftwImpl::OutAlgoDataType; + using InUniquePtr = typename fftwImpl::InUniquePtr; + using OutUniquePtr = typename fftwImpl::OutUniquePtr; + using PlanUniquePtr = typename fftwImpl::PlanUniquePtr; + + std::size_t fftSize{ 0 }; + std::string wisdomPath{ ".gr_fftw_wisdom" }; + int sign{ FFTW_FORWARD }; + unsigned int flags{ FFTW_ESTIMATE }; // FFTW_EXHAUSTIVE, FFTW_MEASURE, FFTW_ESTIMATE + InUniquePtr fftwIn{}; + OutUniquePtr fftwOut{}; + PlanUniquePtr fftwPlan{}; + + FFTw() = default; + FFTw(const FFTw &rhs) = delete; + FFTw(FFTw &&rhs) noexcept = delete; + FFTw & + operator=(const FFTw &rhs) + = delete; + FFTw & + operator=(FFTw &&rhs) noexcept + = delete; + + ~FFTw() { clearFftw(); } + + std::vector + computeFFT(const std::vector &in) noexcept { + if (fftSize != in.size()) { + fftSize = in.size(); + initAll(); + } + std::vector out(getOutputSize()); + computeFFT(in, out); + return out; + } + + void + computeFFT(const std::vector &in, std::vector &out) noexcept { + if (!std::has_single_bit(in.size())) { + throw std::runtime_error(fmt::format("Input data must have 2^N samples, input size: ", in.size())); + } + + if (fftSize != in.size()) { + fftSize = in.size(); + initAll(); + } + + if (out.size() < getOutputSize()) { + throw std::runtime_error(fmt::format("Output vector size ({}) is not enough, at least {} needed. ", out.size(), getOutputSize())); + } + + static_assert(sizeof(InAlgoDataType) == sizeof(T), "Input std::complex and T[2] must have the same size."); + std::memcpy(fftwIn.get(), &(*in.begin()), sizeof(InAlgoDataType) * fftSize); + + fftwImpl::execute(fftwPlan.get()); + + static_assert(sizeof(OutDataType) == sizeof(OutAlgoDataType), "Output std::complex and T[2] must have the same size."); + std::memcpy(out.data(), fftwOut.get(), sizeof(OutDataType) * getOutputSize()); + } + + int + importWisdom() { + // TODO: lock file while importing wisdom? + return fftwImpl::importWisdomFromFilename(wisdomPath); + } + + int + exportWisdom() { + // TODO: lock file while exporting wisdom? + return fftwImpl::exportWisdomToFilename(wisdomPath); + } + + int + importWisdomFromString(const std::string wisdomString) { + return fftwImpl::importWisdomFromString(wisdomString); + } + + std::string + exportWisdomToString() { + return fftwImpl::exportWisdomToString(); + } + + void + forgetWisdom() { + return fftwImpl::forgetWisdom(); + } + +private: + constexpr std::size_t + getOutputSize() { + if constexpr (ComplexType) { + return fftSize; + } else { + return 1 + fftSize / 2; + } + } + + void + initAll() { + clearFftw(); + fftwIn = InUniquePtr(static_cast(fftwImpl::malloc(sizeof(InAlgoDataType) * fftSize))); + fftwOut = OutUniquePtr(static_cast(fftwImpl::malloc(sizeof(OutAlgoDataType) * getOutputSize()))); + + { + std::lock_guard lg{ fftw_plan_mutex }; + importWisdom(); + fftwPlan = PlanUniquePtr(fftwImpl::plan(static_cast(fftSize), fftwIn.get(), fftwOut.get(), sign, flags)); + exportWisdom(); + } + } + + void + clearFftw() { + { + std::lock_guard lg{ fftw_plan_mutex }; + fftwPlan.reset(); + } + fftwIn.reset(); + fftwOut.reset(); + // fftw::cleanup(); // No need for fftw_cleanup -> After calling it, all existing plans become undefined + } +}; + +} // namespace gr::algorithm + +#endif // GRAPH_PROTOTYPE_ALGORITHM_FFTW_HPP diff --git a/test/blocklib/core/fft/window.hpp b/test/blocklib/algorithm/fft/window.hpp similarity index 91% rename from test/blocklib/core/fft/window.hpp rename to test/blocklib/algorithm/fft/window.hpp index 70df9b192..d8831c27d 100644 --- a/test/blocklib/core/fft/window.hpp +++ b/test/blocklib/algorithm/fft/window.hpp @@ -1,24 +1,25 @@ -#ifndef GRAPH_PROTOTYPE_WINDOW_HPP -#define GRAPH_PROTOTYPE_WINDOW_HPP +#ifndef GRAPH_PROTOTYPE_ALGORITHM_WINDOW_HPP +#define GRAPH_PROTOTYPE_ALGORITHM_WINDOW_HPP +#include #include #include #include #include -namespace gr::blocks::fft { +namespace gr::algorithm { -template -concept FloatOrDoubleType = std::is_same_v || std::is_same_v; - -// Implementation of window function (also known as an apodization function or tapering function). -// See Wikipedia for more info: https://en.wikipedia.org/wiki/Window_function +/* + * Implementation of window function (also known as an apodization function or tapering function). + * See Wikipedia for more info: https://en.wikipedia.org/wiki/Window_function + */ enum class WindowFunction : int { None, Rectangular, Hamming, Hann, HannExp, Blackman, Nuttall, BlackmanHarris, BlackmanNuttall, FlatTop, Exponential }; // Only float or double are allowed because FFT supports only float or double precisions. -template +template + requires(std::is_floating_point_v) std::vector createWindowFunction(WindowFunction func, const std::size_t n) { constexpr T pi = std::numbers::pi_v; @@ -110,6 +111,6 @@ createWindowFunction(WindowFunction func, const std::size_t n) { } } } -} // namespace gr::blocks::fft +} // namespace gr::algorithm -#endif // GRAPH_PROTOTYPE_WINDOW_HPP +#endif // GRAPH_PROTOTYPE_ALGORITHM_WINDOW_HPP diff --git a/test/blocklib/core/fft/fft.hpp b/test/blocklib/core/fft/fft.hpp index a9920cfd9..0f2afe196 100644 --- a/test/blocklib/core/fft/fft.hpp +++ b/test/blocklib/core/fft/fft.hpp @@ -1,170 +1,89 @@ #ifndef GRAPH_PROTOTYPE_FFT_HPP #define GRAPH_PROTOTYPE_FFT_HPP -#include "window.hpp" +#include "../../algorithm/fft/fft.hpp" +#include "../../algorithm/fft/fft_common.hpp" +#include "../../algorithm/fft/fft_types.hpp" +#include "../../algorithm/fft/fftw.hpp" +#include "../../algorithm/fft/window.hpp" #include -#include +#include #include #include namespace gr::blocks::fft { using namespace fair::graph; +using gr::algorithm::ComplexType; +using gr::algorithm::FFTInternal; -template -concept ComplexType = std::is_same_v> || std::is_same_v>; - -template -concept DoubleType = std::is_same_v> || std::is_same_v; - -template -concept LongDoubleType = std::is_same_v> || std::is_same_v; - -template -concept FFTSupportedTypes = ComplexType || std::integral || std::floating_point; - -inline static std::mutex fftw_plan_mutex; - -template -struct fft : node> { +template::InDataType>> + requires(ComplexType || std::integral || std::floating_point) +struct FFT : node> { public: - static_assert(not LongDoubleType, "long double is not supported"); - - // clang-format off - template - struct fftw { - using PlanType = fftwf_plan; - using InDataType = std::conditional_t, fftwf_complex, float>; - using OutDataType = fftwf_complex; - using PrecisionType = float; - using InHistoryType = std::conditional_t, std::complex, float>; // history_buffer can not store c-arrays -> use std::complex - using InUniquePtr = std::unique_ptr; - using OutUniquePtr = std::unique_ptr; - using PlanUniquePtr = std::unique_ptr, decltype([](PlanType ptr) { fftwf_destroy_plan(ptr); })>; - - static void execute(const PlanType p) { fftwf_execute(p); } - static void cleanup() { fftwf_cleanup(); } - static void * malloc(std::size_t n) { return fftwf_malloc(n);} - static PlanType plan(int p_n, InDataType *p_in, OutDataType *p_out, int p_sign, unsigned int p_flags) { - if constexpr (std::is_same_v) { - return fftwf_plan_dft_r2c_1d(p_n, p_in, p_out, p_flags); - } else { - return fftwf_plan_dft_1d(p_n, p_in, p_out, p_sign, p_flags); - } - } - static int importWisdomFromFilename(const std::string& path) {return fftwf_import_wisdom_from_filename(path.c_str());} - static int exportWisdomToFilename(const std::string& path) {return fftwf_export_wisdom_to_filename(path.c_str());} - static int importWisdomFromString(const std::string& str) {return fftwf_import_wisdom_from_string(str.c_str());} - static std::string exportWisdomToString() { - char* cstr = fftwf_export_wisdom_to_string(); - if (cstr == nullptr) return ""; - std::string str(cstr); - fftwf_free(cstr); - return str; - } - static void forgetWisdom() {fftwf_forget_wisdom();} - }; - - template - struct fftw { - using PlanType = fftw_plan; - using InDataType = std::conditional_t, fftw_complex, double>; - using OutDataType = fftw_complex; - using PrecisionType = double; - using InHistoryType = std::conditional_t, std::complex, double>; // history_buffer can not store c-arrays -> use std::complex - using InUniquePtr = std::unique_ptr; - using OutUniquePtr = std::unique_ptr; - using PlanUniquePtr = std::unique_ptr, decltype([](PlanType ptr) { fftw_destroy_plan(ptr); })>; - - static void execute(const PlanType p) { fftw_execute(p); } - static void cleanup() { fftw_cleanup(); } - static void * malloc(std::size_t n) { return fftw_malloc(n);} - static PlanType plan(int p_n, InDataType *p_in, OutDataType *p_out, int p_sign, unsigned int p_flags) { - if constexpr (std::is_same_v) { - return fftw_plan_dft_r2c_1d(p_n, p_in, p_out, p_flags); - } else { - return fftw_plan_dft_1d(p_n, p_in, p_out, p_sign, p_flags); - } - } - static int importWisdomFromFilename(const std::string& path) {return fftw_import_wisdom_from_filename(path.c_str());} - static int exportWisdomToFilename(const std::string& path) {return fftw_export_wisdom_to_filename(path.c_str());} - static int importWisdomFromString(const std::string& str) {return fftw_import_wisdom_from_string(str.c_str());} - static std::string exportWisdomToString() { - char* cstr = fftw_export_wisdom_to_string(); - if (cstr == nullptr) return ""; - std::string str(cstr); - fftw_free(cstr); - return str; - } - static void forgetWisdom() {fftw_forget_wisdom();} + using U = gr::algorithm::FFTOutputType; + using InDataType = typename FFTInternal::InDataType; + using OutDataType = typename FFTInternal::OutDataType; + using PrecisionType = typename FFTInternal::PrecisionType; + + PortIn in; + PortOut> out; + + FourierAlgoImpl fftImpl; + Annotated> fftSize{ 1024 }; + Annotated> window{ static_cast(gr::algorithm::WindowFunction::Hann) }; + Annotated, Unit<"Hz">> sampleRate = 1.f; + Annotated signalName = std::string("unknown signal"); + Annotated> signalUnit = std::string("a.u."); + Annotated> signalMin = std::numeric_limits::lowest(); + Annotated> signalMax = std::numeric_limits::max(); + Annotated> outputInDb{ false }; + std::vector windowVector; + history_buffer inputHistory{ fftSize }; + std::vector inData{}; + std::vector outData{}; + std::vector magnitudeSpectrum{}; + std::vector phaseSpectrum{}; + + FFT() { + initAll(); + initWindowFunction(); }; - template struct MakeSigned { using type = signedT;}; - template struct MakeSigned { using type = std::make_signed_t; }; - template struct EvalOutputType { using type1 = evalU; using type2 = typename MakeSigned::type;}; - template struct EvalOutputType { using type1 = typename evalU::value_type; using type2 = evalU;}; - using U = std::conditional_t, typename EvalOutputType::type1, typename EvalOutputType::type2>; - // clang-format on - - using InDataType = typename fftw::InDataType; - using OutDataType = typename fftw::OutDataType; - using PrecisionType = typename fftw::PrecisionType; - using InHistoryType = typename fftw::InHistoryType; - using InUniquePtr = typename fftw::InUniquePtr; - using OutUniquePtr = typename fftw::OutUniquePtr; - using PlanUniquePtr = typename fftw::PlanUniquePtr; + FFT(std::initializer_list> init_parameter) noexcept : node>(init_parameter) {} - PortIn in; - PortOut> out; - - std::size_t fftSize{ 1024 }; - int window{ static_cast(WindowFunction::None) }; - std::vector windowVector; - bool outputInDb{ false }; - std::string wisdomPath{ ".gr_fftw_wisdom" }; - int sign{ FFTW_FORWARD }; - unsigned int flags{ FFTW_ESTIMATE }; // FFTW_EXHAUSTIVE, FFTW_MEASURE, FFTW_ESTIMATE - history_buffer inputHistory{ fftSize }; - InUniquePtr fftwIn{}; - OutUniquePtr fftwOut{}; - PlanUniquePtr fftwPlan{}; - std::vector magnitudeSpectrum{}; - std::vector phaseSpectrum{}; - - fft() { initAll(); } - - fft(const fft &rhs) = delete; - fft(fft &&rhs) noexcept = delete; - fft & - operator=(const fft &rhs) + FFT(const FFT &rhs) = delete; + FFT(FFT &&rhs) noexcept = delete; + FFT & + operator=(const FFT &rhs) = delete; - fft & - operator=(fft &&rhs) noexcept + FFT & + operator=(FFT &&rhs) noexcept = delete; - ~fft() { clearFftw(); } + ~FFT() { clear(); } void settings_changed(const property_map & /*old_settings*/, const property_map &newSettings) noexcept { if (newSettings.contains("fftSize") && fftSize != inputHistory.capacity()) { - inputHistory = history_buffer(fftSize); initAll(); - } else if (newSettings.contains("window")) { // no need to create window twice if fftSize was changed + initWindowFunction(); + } else if (newSettings.contains("window")) { initWindowFunction(); } } constexpr DataSet process_one(T input) noexcept { - if constexpr (std::is_same_v) { + if constexpr (std::is_same_v) { inputHistory.push_back(input); } else { - inputHistory.push_back(static_cast(input)); + inputHistory.push_back(static_cast(input)); } if (inputHistory.size() >= fftSize) { prepareInput(); - computeFft(); + computeFFT(); computeMagnitudeSpectrum(); computePhaseSpectrum(); return createDataset(); @@ -178,161 +97,119 @@ struct fft : node> { DataSet ds{}; ds.timestamp = 0; const std::size_t N{ magnitudeSpectrum.size() }; - - ds.axis_names = { "time", "fft.real", "fft.imag", "magnitude", "phase" }; - ds.axis_units = { "u.a.", "u.a.", "u.a.", "u.a.", "rad" }; - ds.extents = { 4, static_cast(N) }; - ds.layout = fair::graph::layout_right{}; - ds.signal_names = { "fft.real", "fft.imag", "magnitude", "phase" }; - ds.signal_units = { "u.a.", "u.a.", "u.a.", "rad" }; - - ds.signal_values.resize(4 * N); - ds.signal_ranges = { { std::numeric_limits::max(), std::numeric_limits::lowest() }, - { std::numeric_limits::max(), std::numeric_limits::lowest() }, - { std::numeric_limits::max(), std::numeric_limits::lowest() }, - { std::numeric_limits::max(), std::numeric_limits::lowest() } }; - for (std::size_t i = 0; i < N; i++) { - if constexpr (std::is_same_v) { - ds.signal_values[i] = fftwOut[i][0]; - ds.signal_values[i + N] = fftwOut[i][1]; + const std::size_t dim = 5; + + ds.axis_names = { "Frequency", "Re(FFT)", "Im(FFT)", "Magnitude", "Phase" }; + ds.axis_units = { "Hz", signalUnit, fmt::format("i{}", signalUnit), fmt::format("{}/√Hz", signalUnit), "rad" }; + ds.extents = { dim, static_cast(N) }; + ds.layout = fair::graph::layout_right{}; + ds.signal_names = { signalName, fmt::format("Re(FFT({}))", signalName), fmt::format("Im(FFT({}))", signalName), fmt::format("Magnitude({})", signalName), fmt::format("Phase({})", signalName) }; + ds.signal_units = { "Hz", signalUnit, fmt::format("i{}", signalUnit), fmt::format("{}/√Hz", signalUnit), "rad" }; + + ds.signal_values.resize(dim * N); + std::ranges::transform(std::views::iota(std::size_t(0), N), std::ranges::begin(ds.signal_values), [s = sampleRate, f = fftSize, N](const auto i) { + auto const freq = static_cast(s) / static_cast(f); + auto const freqi = static_cast(i) * freq; + if constexpr (ComplexType) { + return freqi - static_cast(N / 2) * freq; } else { - ds.signal_values[i] = static_cast(fftwOut[i][0]); - ds.signal_values[i + N] = static_cast(fftwOut[i][1]); + return freqi; } - ds.signal_ranges[0][0] = std::min(ds.signal_ranges[0][0], ds.signal_values[i]); - ds.signal_ranges[0][1] = std::max(ds.signal_ranges[0][1], ds.signal_values[i]); - ds.signal_ranges[1][0] = std::min(ds.signal_ranges[1][0], ds.signal_values[i + N]); - ds.signal_ranges[1][1] = std::max(ds.signal_ranges[1][1], ds.signal_values[i + N]); - - ds.signal_values[i + 2 * N] = magnitudeSpectrum[i]; - ds.signal_ranges[2][0] = std::min(ds.signal_ranges[2][0], magnitudeSpectrum[i]); - ds.signal_ranges[2][1] = std::max(ds.signal_ranges[2][1], magnitudeSpectrum[i]); - - ds.signal_values[i + 3 * N] = phaseSpectrum[i]; - ds.signal_ranges[3][0] = std::min(ds.signal_ranges[3][0], phaseSpectrum[i]); - ds.signal_ranges[3][1] = std::max(ds.signal_ranges[3][1], phaseSpectrum[i]); + }); + std::transform(outData.begin(), outData.end(), std::next(ds.signal_values.begin(), N), [](const auto &c) { return c.real(); }); + std::transform(outData.begin(), outData.end(), std::next(ds.signal_values.begin(), 2U * N), [](const auto &c) { return c.imag(); }); + std::copy_n(magnitudeSpectrum.begin(), N, std::next(ds.signal_values.begin(), 3U * N)); + std::copy_n(phaseSpectrum.begin(), N, std::next(ds.signal_values.begin(), 4U * N)); + + ds.signal_ranges.resize(dim); + for (std::size_t i = 0; i < dim; i++) { + const auto mm = std::minmax_element(std::next(ds.signal_values.begin(), static_cast(i * N)), std::next(ds.signal_values.begin(), static_cast((i + 1U) * N))); + ds.signal_ranges[i] = { *mm.first, *mm.second }; } ds.signal_errors = {}; - ds.meta_information = { - { { "fftSize", fftSize }, { "window", window }, { "outputInDb", outputInDb }, { "numerator", this->numerator }, { "denominator", this->denominator }, { "stride", this->stride } } - }; + ds.meta_information = { { { "sampleRate", sampleRate }, + { "signalName", signalName }, + { "signalUnit", signalUnit }, + { "signalMin", signalMin }, + { "signalMax", signalMax }, + { "fftSize", fftSize }, + { "window", window }, + { "outputInDb", outputInDb }, + { "numerator", this->numerator }, + { "denominator", this->denominator }, + { "stride", this->stride } } }; return ds; } void prepareInput() { - static_assert(sizeof(InDataType) == sizeof(InHistoryType), "std::complex and T[2] must have the same size"); - std::memcpy(fftwIn.get(), &(*inputHistory.begin()), sizeof(InDataType) * fftSize); + std::copy_n(inputHistory.begin(), fftSize, inData.begin()); // apply window function if needed - if (window != static_cast(WindowFunction::None)) { - if (fftSize != windowVector.size()) throw std::runtime_error(fmt::format("fftSize({}) and windowVector.size({}) are not equal.", fftSize, windowVector.size())); + if (window != static_cast(gr::algorithm::WindowFunction::None)) { + if (fftSize != windowVector.size()) { + throw std::runtime_error(fmt::format("fftSize({}) and windowVector.size({}) are not equal.", fftSize, windowVector.size())); + } for (std::size_t i = 0; i < fftSize; i++) { if constexpr (ComplexType) { - fftwIn[i][0] *= windowVector[i]; - fftwIn[i][1] *= windowVector[i]; + inData[i].real(inData[i].real() * windowVector[i]); + inData[i].imag(inData[i].imag() * windowVector[i]); } else { - fftwIn[i] *= windowVector[i]; + inData[i] *= windowVector[i]; } } } } - void - computeFft() { - fftw::execute(fftwPlan.get()); + inline void + computeFFT() { + // fftImpl.computeFFT(inData, outData); + outData = fftImpl.computeFFT(inData); } - void + inline void computeMagnitudeSpectrum() { - for (std::size_t i = 0; i < magnitudeSpectrum.size(); i++) { - const PrecisionType mag{ std::hypot(fftwOut[i][0], fftwOut[i][1]) * static_cast(2.0) / static_cast(fftSize) }; - magnitudeSpectrum[i] = static_cast(outputInDb ? static_cast(20.) * std::log10(std::abs(mag)) : mag); - } + gr::algorithm::computeMagnitudeSpectrum(outData, magnitudeSpectrum, fftSize, outputInDb); } - void + inline void computePhaseSpectrum() { - for (std::size_t i = 0; i < phaseSpectrum.size(); i++) { - const auto phase{ std::atan2(fftwOut[i][1], fftwOut[i][0]) }; - phaseSpectrum[i] = outputInDb ? static_cast(20.) * static_cast(std::log10(std::abs(phase))) : static_cast(phase); - } + gr::algorithm::computePhaseSpectrum(outData, phaseSpectrum); } void initAll() { - clearFftw(); - fftwIn = InUniquePtr(static_cast(fftw::malloc(sizeof(InDataType) * fftSize))); + clear(); + inputHistory = history_buffer(fftSize); + inData.resize(fftSize); if constexpr (ComplexType) { - fftwOut = OutUniquePtr(static_cast(fftw::malloc(sizeof(OutDataType) * fftSize))); + outData.resize(fftSize); magnitudeSpectrum.resize(fftSize); phaseSpectrum.resize(fftSize); } else { - fftwOut = OutUniquePtr(static_cast(fftw::malloc(sizeof(OutDataType) * (1 + fftSize / 2)))); + outData.resize(1 + fftSize / 2); magnitudeSpectrum.resize(fftSize / 2); phaseSpectrum.resize(fftSize / 2); } - { - std::lock_guard lg{ fftw_plan_mutex }; - importWisdom(); - fftwPlan = PlanUniquePtr(fftw::plan(static_cast(fftSize), fftwIn.get(), fftwOut.get(), sign, flags)); - exportWisdom(); - } - - initWindowFunction(); } - void + inline void initWindowFunction() { - windowVector = createWindowFunction(static_cast(window), fftSize); + windowVector = createWindowFunction(static_cast(window.value), fftSize); } void - clearFftw() { - { - std::lock_guard lg{ fftw_plan_mutex }; - fftwPlan.reset(); - } - - fftwIn.reset(); - fftwOut.reset(); - - // fftw::cleanup(); // No need for fftw_cleanup -> After calling it, all existing plans become undefined + clear() { + inData.clear(); + outData.clear(); magnitudeSpectrum.clear(); phaseSpectrum.clear(); } - - int - importWisdom() { - // TODO: lock file while importing wisdom? - return fftw::importWisdomFromFilename(wisdomPath); - } - - int - exportWisdom() { - // TODO: lock file while exporting wisdom? - return fftw::exportWisdomToFilename(wisdomPath); - } - - int - importWisdomFromString(const std::string wisdomString) { - return fftw::importWisdomFromString(wisdomString); - } - - std::string - exportWisdomToString() { - return fftw::exportWisdomToString(); - } - - void - forgetWisdom() { - return fftw::forgetWisdom(); - } }; } // namespace gr::blocks::fft -ENABLE_REFLECTION_FOR_TEMPLATE_FULL((typename T), (gr::blocks::fft::fft), in, out, fftSize, outputInDb, window); +ENABLE_REFLECTION_FOR_TEMPLATE_FULL((typename T, typename TImpl), (gr::blocks::fft::FFT), in, out, fftSize, sampleRate, signalName, signalUnit, signalMin, signalMax, outputInDb, window); #endif // GRAPH_PROTOTYPE_FFT_HPP diff --git a/test/qa_fft.cpp b/test/qa_fft.cpp index 87b978da9..fa2bdf31e 100644 --- a/test/qa_fft.cpp +++ b/test/qa_fft.cpp @@ -6,6 +6,8 @@ template<> auto boost::ut::cfg = boost::ut::runner>{}; #endif +#include "blocklib/algorithm/fft/fft.hpp" +#include "blocklib/algorithm/fft/fftw.hpp" #include "blocklib/core/fft/fft.hpp" #include #include @@ -63,25 +65,86 @@ equalVectors(const std::vector &v1, const std::vector &v2, double toleranc template void -testFftwTypes() { +testFFTwTypes() { using namespace boost::ut; - gr::blocks::fft::fft fft1; + gr::algorithm::FFTw fft1; expect(std::is_same_v, inT>) << ""; expect(std::is_same_v, outT>) << ""; expect(std::is_same_v) << ""; } +template +void +equalDataset(const gr::blocks::fft::FFT &fft1, const fair::graph::DataSet &ds1, float sampleRate) { + using namespace boost::ut; + using namespace boost::ut::reflection; + + const U tolerance = U(0.0001); + + const auto N = fft1.magnitudeSpectrum.size(); + if (N == fft1.fftSize) { // complex input + // fmt::println("ds: {} exp:{}", ds1.signal_values[0], -(static_cast(N) / 2.) * sampleRate); + // expect(approx(ds1.signal_values[0], -(static_cast(N) / 2.) * sampleRate, tolerance)) << fmt::format("<{}> equal DataSet frequency[0]", type_name()); + + } else { // real input + }; + bool isEqualFFTOut = true; + const auto NSize = static_cast(N); + for (std::size_t i = 0; i < NSize; i++) { + if (ds1.signal_values[i + NSize] != static_cast(fft1.outData[i].real()) || ds1.signal_values[i + 2U * NSize] != static_cast(fft1.outData[i].imag())) { + isEqualFFTOut = false; + break; + } + } + expect(eq(isEqualFFTOut, true)) << fmt::format("<{}> equal DataSet FFT output", type_name()); + expect(equalVectors(std::vector(ds1.signal_values.begin() + static_cast(3U * N), ds1.signal_values.begin() + static_cast(4U * N)), fft1.magnitudeSpectrum)) + << fmt::format("<{}> equal DataSet magnitude", type_name()); + expect(equalVectors(std::vector(ds1.signal_values.begin() + static_cast(4U * N), ds1.signal_values.begin() + static_cast(5U * N)), fft1.phaseSpectrum)) + << fmt::format("<{}> equal DataSet phase", type_name()); + + for (std::size_t i = 0; i < 5; i++) { + const auto mm = std::minmax_element(std::next(ds1.signal_values.begin(), static_cast(i * N)), std::next(ds1.signal_values.begin(), static_cast((i + 1U) * N))); + if constexpr (std::is_integral_v) { + expect(eq(*mm.first, ds1.signal_ranges[i][0])); + expect(eq(*mm.second, ds1.signal_ranges[i][1])); + } else { + constexpr U tolerance{ static_cast(0.000001) }; + expect(approx(*mm.first, ds1.signal_ranges[i][0], tolerance)); + expect(approx(*mm.second, ds1.signal_ranges[i][1], tolerance)); + } + } +} + +template +using FFTwAlgo = gr::algorithm::FFTw::InDataType>; + +template +using FFTAlgo = gr::algorithm::FFT::InDataType>; + +template typename A> +struct TypePair { + using DataType = T; + using AlgoType = A; +}; + const boost::ut::suite fftTests = [] { using namespace boost::ut; using namespace gr::blocks::fft; using namespace boost::ut::reflection; - std::tuple, std::complex> complexTypesToTest{}; - std::tuple, std::complex, float, double, int> typesToTest{}; - std::tuple floatingTypesToTest{}; + + std::tuple, FFTwAlgo>, TypePair, FFTAlgo>, TypePair, FFTwAlgo>, TypePair, FFTAlgo>> + complexTypesWithAlgoToTest{}; + std::tuple, FFTwAlgo>, TypePair, FFTAlgo>, TypePair, FFTwAlgo>, TypePair, FFTAlgo>, + TypePair, TypePair, TypePair, TypePair, TypePair, TypePair> + typesWithAlgoToTest{}; + + std::tuple floatingTypesToTest{}; "FFT sin tests"_test = []() { - fft fft1{}; - constexpr double tolerance{ 1.e-6 }; + using DataType = T::DataType; + using AlgoType = T::AlgoType; + FFT fft1{}; + constexpr double tolerance{ 1.e-5 }; struct TestParams { std::size_t N{ 1024 }; // must be power of 2 double sampleRate{ 128. }; // must be power of 2 (only for the unit test for easy comparison with true result) @@ -97,13 +160,14 @@ const boost::ut::suite fftTests = [] { std::ignore = fft1.settings().set({ { "fftSize", t.N } }); std::ignore = fft1.settings().set({ { "outputInDb", t.outputInDb } }); + std::ignore = fft1.settings().set({ { "window", static_cast(gr::algorithm::WindowFunction::None) } }); std::ignore = fft1.settings().apply_staged_parameters(); - const auto signal{ generateSinSample(t.N, t.sampleRate, t.frequency, t.amplitude) }; + const auto signal{ generateSinSample(t.N, t.sampleRate, t.frequency, t.amplitude) }; for (const auto s : signal) { - fft1.inputHistory.push_back(static_cast::InHistoryType>(s)); + fft1.inputHistory.push_back(static_cast::InDataType>(s)); } fft1.prepareInput(); - fft1.computeFft(); + fft1.computeFFT(); fft1.computeMagnitudeSpectrum(); const auto peakIndex{ @@ -114,37 +178,38 @@ const boost::ut::suite fftTests = [] { const auto peakFrequency{ static_cast(peakIndex) * t.sampleRate / static_cast(t.N) }; const auto expectedAmplitude = t.outputInDb ? 20. * log10(std::abs(t.amplitude)) : t.amplitude; - if constexpr (!std::is_same_v) { - expect(approx(static_cast(peakAmplitude), expectedAmplitude, tolerance)) << fmt::format("<{}> equal amplitude", type_name()); - expect(approx(peakFrequency, t.frequency, tolerance)) << fmt::format("<{}> equal frequency", type_name()); + if constexpr (!std::is_same_v) { + expect(approx(static_cast(peakAmplitude), expectedAmplitude, tolerance)) << fmt::format("{} equal amplitude", type_name()); + expect(approx(peakFrequency, t.frequency, tolerance)) << fmt::format("{} equal frequency", type_name()); } } - } | typesToTest; + } | typesWithAlgoToTest; "FFT pattern tests"_test = []() { - fft fft1{}; - constexpr double tolerance{ 1.e-6 }; - constexpr std::size_t N{ 16 }; - std::ignore = fft1.settings().set({ { "fftSize", N } }); + using DataType = T::DataType; + using AlgoType = T::AlgoType; + constexpr double tolerance{ 1.e-5 }; + constexpr std::size_t N{ 16 }; + FFT fft1({ { "fftSize", N }, { "window", static_cast(gr::algorithm::WindowFunction::None) } }); std::ignore = fft1.settings().apply_staged_parameters(); - std::vector signal(N); + std::vector signal(N); static_assert(N == 16, "expected values are calculated for N == 16"); std::size_t expectedPeakIndex{ 0 }; - T expectedFft0{ 0., 0. }; + DataType expectedFft0{ 0., 0. }; double expectedPeakAmplitude{ 0. }; for (std::size_t iT = 0; iT < 5; iT++) { if (iT == 0) { - std::fill(signal.begin(), signal.end(), T(0., 0.)); + std::fill(signal.begin(), signal.end(), DataType(0., 0.)); expectedFft0 = { 0., 0. }; expectedPeakAmplitude = 0.; } else if (iT == 1) { - std::fill(signal.begin(), signal.end(), T(1., 0.)); + std::fill(signal.begin(), signal.end(), DataType(1., 0.)); expectedFft0 = { 16., 0. }; expectedPeakAmplitude = 2.; } else if (iT == 2) { - std::fill(signal.begin(), signal.end(), T(1., 1.)); + std::fill(signal.begin(), signal.end(), DataType(1., 1.)); expectedFft0 = { 16., 16. }; expectedPeakAmplitude = std::sqrt(8.); } else if (iT == 3) { @@ -153,14 +218,14 @@ const boost::ut::suite fftTests = [] { expectedPeakAmplitude = 17.; } else if (iT == 4) { int i = 0; - std::generate(signal.begin(), signal.end(), [&i] { return T(static_cast(i++ % 2), 0.); }); + std::generate(signal.begin(), signal.end(), [&i] { return DataType(static_cast(i++ % 2), 0.); }); expectedFft0 = { 8., 0. }; expectedPeakAmplitude = 1.; } fft1.inputHistory.push_back_bulk(signal.begin(), signal.end()); fft1.prepareInput(); - fft1.computeFft(); + fft1.computeFFT(); fft1.computeMagnitudeSpectrum(); const auto peakIndex{ static_cast(std::distance(fft1.magnitudeSpectrum.begin(), std::max_element(fft1.magnitudeSpectrum.begin(), fft1.magnitudeSpectrum.end()))) }; @@ -168,68 +233,55 @@ const boost::ut::suite fftTests = [] { expect(eq(peakIndex, expectedPeakIndex)) << fmt::format("<{}> equal peak index", type_name()); expect(approx(static_cast(peakAmplitude), expectedPeakAmplitude, tolerance)) << fmt::format("<{}> equal amplitude", type_name()); - expect(approx(static_cast(fft1.fftwOut[0][0]), static_cast(expectedFft0.real()), tolerance)) << fmt::format("<{}> equal fft[0].real()", type_name()); - expect(approx(static_cast(fft1.fftwOut[0][1]), static_cast(expectedFft0.imag()), tolerance)) << fmt::format("<{}> equal fft[0].imag()", type_name()); + expect(approx(static_cast(fft1.outData[0].real()), static_cast(expectedFft0.real()), tolerance)) << fmt::format("<{}> equal fft[0].real()", type_name()); + expect(approx(static_cast(fft1.outData[0].imag()), static_cast(expectedFft0.imag()), tolerance)) << fmt::format("<{}> equal fft[0].imag()", type_name()); } - } | complexTypesToTest; + } | complexTypesWithAlgoToTest; "FFT process_one tests"_test = []() { - fft fft1{}; - constexpr std::size_t N{ 16 }; - std::ignore = fft1.settings().set({ { "fftSize", N } }); - std::ignore = fft1.settings().apply_staged_parameters(); - using DatasetType = typename fft::U; - using InHistoryType = typename fft::InHistoryType; - - std::vector signal(N); + using DataType = T::DataType; + using AlgoType = T::AlgoType; + using DatasetType = typename FFT::U; + using InDataType = typename FFT::InDataType; + + constexpr std::size_t N{ 16 }; + constexpr float sampleRate{ 1.f }; + FFT fft1({ { "fftSize", N }, { "sampleRate", sampleRate } }); + std::ignore = fft1.settings().apply_staged_parameters(); + + std::vector signal(N); std::iota(signal.begin(), signal.end(), 1); DataSet ds1{}; - for (std::size_t i = 0; i < N; i++) ds1 = fft1.process_one(signal[i]); - expect(equalVectors(std::vector(fft1.inputHistory.begin(), fft1.inputHistory.end()), signal)) << fmt::format("<{}> equal history buffer", type_name()); - const auto N2 = fft1.magnitudeSpectrum.size(); - expect(equalVectors(std::vector(ds1.signal_values.begin() + static_cast(2U * N2), ds1.signal_values.begin() + static_cast(3U * N2)), - fft1.magnitudeSpectrum)) - << fmt::format("<{}> equal DataSet magnitude", type_name()); - - for (std::size_t i = 0; i < 4; i++) { - const auto mm = std::minmax_element(std::next(ds1.signal_values.begin(), static_cast(i * N2)), - std::next(ds1.signal_values.begin(), static_cast((i + 1U) * N2))); - if constexpr (std::is_integral_v) { - expect(eq(*mm.first, ds1.signal_ranges[i][0])); - expect(eq(*mm.second, ds1.signal_ranges[i][1])); - } else { - constexpr DatasetType tolerance{ static_cast(0.000001) }; - expect(approx(*mm.first, ds1.signal_ranges[i][0], tolerance)); - expect(approx(*mm.second, ds1.signal_ranges[i][1], tolerance)); - } + for (const auto s : signal) { + ds1 = fft1.process_one(s); } + expect(equalVectors(std::vector(fft1.inputHistory.begin(), fft1.inputHistory.end()), signal)) << fmt::format("<{}> equal history buffer", type_name()); + equalDataset(fft1, ds1, sampleRate); std::iota(signal.begin(), signal.end(), N + 1); - for (std::size_t i = 0; i < N; i++) ds1 = fft1.process_one(signal[i]); - expect(equalVectors(std::vector(fft1.inputHistory.begin(), fft1.inputHistory.end()), signal)) << fmt::format("<{}> equal history buffer", type_name()); - expect(equalVectors(std::vector(ds1.signal_values.begin() + static_cast(2U * N2), ds1.signal_values.begin() + static_cast(3U * N2)), - fft1.magnitudeSpectrum)) - << fmt::format("<{}> equal DataSet magnitude", type_name()); - } | typesToTest; + for (std::size_t i = 0; i < N; i++) { + ds1 = fft1.process_one(signal[i]); + } + expect(equalVectors(std::vector(fft1.inputHistory.begin(), fft1.inputHistory.end()), signal)) << fmt::format("<{}> equal history buffer", type_name()); + equalDataset(fft1, ds1, sampleRate); + } | typesWithAlgoToTest; "FFT types tests"_test = [] { - expect(std::is_same_v>::U, float>) << "output type must be float"; - expect(std::is_same_v>::U, double>) << "output type must be double"; - expect(std::is_same_v::U, float>) << "output type must be float"; - expect(std::is_same_v::U, double>) << "output type must be double"; - expect(std::is_same_v::U, int>) << "output type must be int"; - expect(std::is_same_v::U, int>) << "output type must be int"; - expect(std::is_same_v::U, int64_t>) << "output type must be int64_t"; - expect(std::is_same_v::U, int64_t>) << "output type must be int64_t"; + expect(std::is_same_v>::U, float>) << "output type must be float"; + expect(std::is_same_v>::U, double>) << "output type must be double"; + expect(std::is_same_v::U, float>) << "output type must be float"; + expect(std::is_same_v::U, double>) << "output type must be double"; + expect(std::is_same_v::U, int>) << "output type must be int"; + expect(std::is_same_v::U, int>) << "output type must be int"; + expect(std::is_same_v::U, int64_t>) << "output type must be int64_t"; + expect(std::is_same_v::U, int64_t>) << "output type must be int64_t"; }; "FFT fftw types tests"_test = [] { - testFftwTypes, fftwf_complex, fftwf_complex, fftwf_plan>(); - testFftwTypes, fftw_complex, fftw_complex, fftw_plan>(); - testFftwTypes(); - testFftwTypes(); - testFftwTypes(); - testFftwTypes(); + testFFTwTypes, fftwf_complex, fftwf_complex, fftwf_plan>(); + testFFTwTypes, fftw_complex, fftw_complex, fftw_plan>(); + testFFTwTypes(); + testFFTwTypes(); }; "FFT flow graph example"_test = [] { @@ -240,7 +292,7 @@ const boost::ut::suite fftTests = [] { fg::graph flow1; auto &source1 = flow1.make_node>(); - auto &fft1 = flow1.make_node>({ { "fftSize", 16 } }); + auto &fft1 = flow1.make_node>({ { "fftSize", 16 } }); std::ignore = flow1.connect<"out">(source1).to<"in">(fft1); auto sched1 = Scheduler(std::move(flow1), threadPool); @@ -248,7 +300,7 @@ const boost::ut::suite fftTests = [] { for (int i = 0; i < 2; i++) { fg::graph flow2; auto &source2 = flow2.make_node>(); - auto &fft2 = flow2.make_node>({ { "fftSize", 16 } }); + auto &fft2 = flow2.make_node>({ { "fftSize", 16 } }); std::ignore = flow2.connect<"out">(source2).to<"in">(fft2); auto sched2 = Scheduler(std::move(flow2), threadPool); sched2.run_and_wait(); @@ -260,12 +312,15 @@ const boost::ut::suite fftTests = [] { }; "FFT window function tests"_test = []() { - fft fft1{}; - using PrecisionType = typename fft::PrecisionType; - using InHistoryType = typename fft::InHistoryType; - constexpr PrecisionType tolerance{ PrecisionType(0.00001) }; - constexpr std::size_t N{ 8 }; - + using DataType = T::DataType; + using AlgoType = T::AlgoType; + FFT fft1{}; + + using PrecisionType = typename FFT::PrecisionType; + using InDataType = typename FFT::InDataType; + constexpr PrecisionType tolerance{ PrecisionType(0.00001) }; + constexpr std::size_t N{ 8 }; + using gr::algorithm::WindowFunction; std::vector testCases = { WindowFunction::None, WindowFunction::Rectangular, WindowFunction::Hamming, WindowFunction::Hann, WindowFunction::HannExp, WindowFunction::Blackman, WindowFunction::Nuttall, WindowFunction::BlackmanHarris, WindowFunction::BlackmanNuttall, @@ -276,17 +331,17 @@ const boost::ut::suite fftTests = [] { std::ignore = fft1.settings().set({ { "window", static_cast(t) } }); std::ignore = fft1.settings().apply_staged_parameters(); - std::vector signal(N); - if constexpr (ComplexType) { - typename T::value_type i = 0.; + std::vector signal(N); + if constexpr (ComplexType) { + typename DataType::value_type i = 0.; std::generate(signal.begin(), signal.end(), [&i] { - i = i + static_cast(1.); - return T(i, i); + i = i + static_cast(1.); + return DataType(i, i); }); } else { std::iota(signal.begin(), signal.end(), 1.); } - for (const auto s : signal) fft1.inputHistory.push_back(static_cast(s)); + for (const auto s : signal) fft1.inputHistory.push_back(static_cast(s)); fft1.prepareInput(); expect(eq(fft1.fftSize, N)) << fmt::format("<{}> equal fft size", type_name()); @@ -295,20 +350,22 @@ const boost::ut::suite fftTests = [] { std::vector windowFunc = createWindowFunction(t, N); for (std::size_t i = 0; i < N; i++) { - if constexpr (ComplexType) { - const typename T::value_type expValue = (t == WindowFunction::None) ? signal[i].real() : signal[i].real() * static_cast(windowFunc[i]); - expect(approx(fft1.fftwIn.get()[i][0], expValue, tolerance)) << fmt::format("<{}> equal fftwIn complex.real", type_name()); - expect(approx(fft1.fftwIn.get()[i][1], expValue, tolerance)) << fmt::format("<{}> equal fftwIn complex.imag", type_name()); + if constexpr (ComplexType) { + const typename DataType::value_type expValue = (t == WindowFunction::None) ? signal[i].real() : signal[i].real() * static_cast(windowFunc[i]); + expect(approx(fft1.inData[i].real(), expValue, tolerance)) << fmt::format("<{}> equal fftwIn complex.real", type_name()); + expect(approx(fft1.inData[i].imag(), expValue, tolerance)) << fmt::format("<{}> equal fftwIn complex.imag", type_name()); } else { const PrecisionType expValue = (t == WindowFunction::None) ? static_cast(signal[i]) : static_cast(signal[i]) * static_cast(windowFunc[i]); - expect(approx(fft1.fftwIn.get()[i], expValue, tolerance)) << fmt::format("<{}> equal fftwIn", type_name()); + expect(approx(fft1.inData[i], expValue, tolerance)) << fmt::format("<{}> equal fftwIn", type_name()); } } } - } | typesToTest; + } | typesWithAlgoToTest; "FFT window tests"_test = []() { + using gr::algorithm::WindowFunction; + // Expected value for size 8 std::vector None8{}; std::vector Rectangular8{ 1., 1., 1., 1., 1., 1., 1., 1. }; @@ -354,14 +411,14 @@ const boost::ut::suite fftTests = [] { expect(eq(createWindowFunction(WindowFunction::Nuttall, 0).size(), 0u)) << fmt::format("<{}> zero size Nuttall[8] vectors", type_name()); } | floatingTypesToTest; - "FFT wisdom import/export tests"_test = []() { - fft fft1{}; + "FFTW wisdom import/export tests"_test = []() { + gr::algorithm::FFTw fftw1{}; - std::string wisdomString1 = fft1.exportWisdomToString(); - fft1.forgetWisdom(); - int importOk = fft1.importWisdomFromString(wisdomString1); + std::string wisdomString1 = fftw1.exportWisdomToString(); + fftw1.forgetWisdom(); + int importOk = fftw1.importWisdomFromString(wisdomString1); expect(eq(importOk, 1)) << "Wisdom import from string."; - std::string wisdomString2 = fft1.exportWisdomToString(); + std::string wisdomString2 = fftw1.exportWisdomToString(); // lines are not always at the same order thus it is hard to compare // expect(eq(wisdomString1, wisdomString2)) << "Wisdom strings are the same.";