From bf360a2c5fddac35d356b1de2d415e5e203345c9 Mon Sep 17 00:00:00 2001 From: drslebedev Date: Thu, 12 Oct 2023 09:11:19 +0200 Subject: [PATCH] FFT block and algorithms refactoring: * Change signature of the FFT algorithm classs. Add TInput and TOutput template parameters * Precision of the calculation is defined by the TOutput precision * Add options to calculate phase in deg and unwrapped * Clean up FFT unittests --- .../gnuradio-4.0/algorithm/fourier/fft.hpp | 49 +++-- .../algorithm/fourier/fft_common.hpp | 107 +++++++-- .../algorithm/fourier/fft_types.hpp | 28 --- .../gnuradio-4.0/algorithm/fourier/fftw.hpp | 104 +++++---- algorithm/test/qa_algorithm_fourier.cpp | 172 ++++++++++++++- .../include/gnuradio-4.0/fourier/fft.hpp | 69 +++--- blocks/fourier/test/qa_fourier.cpp | 206 ++++-------------- core/benchmarks/bm_fft.cpp | 28 ++- meta/include/gnuradio-4.0/meta/utils.hpp | 4 + 9 files changed, 447 insertions(+), 320 deletions(-) delete mode 100644 algorithm/include/gnuradio-4.0/algorithm/fourier/fft_types.hpp diff --git a/algorithm/include/gnuradio-4.0/algorithm/fourier/fft.hpp b/algorithm/include/gnuradio-4.0/algorithm/fourier/fft.hpp index d15575317..823041533 100644 --- a/algorithm/include/gnuradio-4.0/algorithm/fourier/fft.hpp +++ b/algorithm/include/gnuradio-4.0/algorithm/fourier/fft.hpp @@ -3,19 +3,17 @@ #include -#include "fft_types.hpp" #include "window.hpp" namespace gr::algorithm { -template - requires(ComplexType || std::floating_point) +template, TInput, std::complex>> + requires((gr::meta::complex_like || std::floating_point) && (gr::meta::complex_like) ) struct FFT { - using PrecisionType = FFTAlgoPrecision::type; - using OutDataType = std::conditional_t, T, std::complex>; + using Precision = TOutput::value_type; - std::vector twiddleFactors{}; - std::size_t fftSize{ 0 }; + std::vector twiddleFactors{}; + std::size_t fftSize{ 0 }; FFT() = default; FFT(const FFT &rhs) = delete; @@ -34,15 +32,16 @@ struct FFT { precomputeTwiddleFactors(); } - std::vector - computeFFT(const std::vector &in) { - std::vector out(in.size()); - computeFFT(in, out); - return out; - } + auto + compute(const std::ranges::input_range auto &in, std::ranges::output_range auto &&out) { + if constexpr (requires(std::size_t n) { out.resize(n); }) { + if (out.size() != in.size()) { + out.resize(in.size()); + } + } else { + static_assert(std::tuple_size_v == std::tuple_size_v, "Size mismatch for fixed-size container."); + } - void - computeFFT(const std::vector &in, std::vector &out) { if (!std::has_single_bit(in.size())) { throw std::invalid_argument(fmt::format("Input data must have 2^N samples, input size: ", in.size())); } @@ -52,9 +51,10 @@ struct FFT { } // For the moment no optimization for real data inputs, just create complex with zero imaginary value. - if constexpr (!ComplexType) { - std::ranges::transform(in.begin(), in.end(), out.begin(), [](const auto c) { return OutDataType(c, 0.); }); + if constexpr (!gr::meta::complex_like) { + std::ranges::transform(in.begin(), in.end(), out.begin(), [](const auto c) { return TOutput(c, 0.); }); } else { + // precision is defined by output type, if cast is needed, let `std::copy` do the job std::ranges::copy(in.begin(), in.end(), out.begin()); } @@ -77,11 +77,18 @@ struct FFT { } } } + + return out; + } + + auto + compute(const std::ranges::input_range auto &in) { + return compute(in, std::vector(in.size())); } private: void - bitReversalPermutation(std::vector &vec) const noexcept { + 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); @@ -92,12 +99,12 @@ struct FFT { void precomputeTwiddleFactors() { twiddleFactors.clear(); - const auto minus2Pi = PrecisionType(-2. * std::numbers::pi); + const auto minus2Pi = Precision(-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))) }; + const TOutput w{ std::exp(TOutput(0., minus2Pi / static_cast(s))) }; for (std::size_t k = 0; k < fftSize; k += s) { - OutDataType wk{ 1., 0. }; + TOutput wk{ 1., 0. }; for (std::size_t j = 0; j < m; j++) { twiddleFactors.push_back(wk); wk *= w; diff --git a/algorithm/include/gnuradio-4.0/algorithm/fourier/fft_common.hpp b/algorithm/include/gnuradio-4.0/algorithm/fourier/fft_common.hpp index ce52cdc8b..23c7e4c23 100644 --- a/algorithm/include/gnuradio-4.0/algorithm/fourier/fft_common.hpp +++ b/algorithm/include/gnuradio-4.0/algorithm/fourier/fft_common.hpp @@ -3,36 +3,109 @@ #include #include +#include #include #include +#include -#include "fft_types.hpp" +namespace gr::algorithm::fft { -namespace gr::algorithm { +struct ConfigMagnitude { + bool computeHalfSpectrum = false; + bool outputInDb = false; +}; -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())); +template TContainerOut = std::vector, + typename T = TContainerIn::value_type> + requires(std::is_same_v> || std::is_same_v>) +auto +computeMagnitudeSpectrum(const TContainerIn &fftIn, TContainerOut &&magOut = {}, ConfigMagnitude config = {}) { + std::size_t magSize = config.computeHalfSpectrum ? fftIn.size() : (fftIn.size() / 2); + if constexpr (requires(std::size_t n) { magOut.resize(n); }) { + if (magOut.size() != magSize) { + magOut.resize(magSize); + } + } else { + static_assert(std::tuple_size_v == std::tuple_size_v, "Size mismatch for fixed-size container."); } + using PrecisionType = typename T::value_type; - std::transform(fftOut.begin(), std::next(fftOut.begin(), static_cast(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); + std::transform(fftIn.begin(), std::next(fftIn.begin(), static_cast(magSize)), magOut.begin(), [fftSize = fftIn.size(), outputInDb = config.outputInDb](const auto &c) { + const auto mag{ std::hypot(c.real(), c.imag()) * PrecisionType(2.) / static_cast(fftSize) }; + return outputInDb ? PrecisionType(20.) * std::log10(std::abs(mag)) : mag; }); + + return magOut; +} + +template + requires(std::is_same_v> || std::is_same_v>) +auto +computeMagnitudeSpectrum(const TContainerIn &fftIn, ConfigMagnitude config) { + return computeMagnitudeSpectrum(fftIn, {}, config); } -template +struct ConfigPhase { + bool computeHalfSpectrum = false; + bool outputInDeg = false; + bool unwrapPhase = false; +}; + +template + requires(std::floating_point) 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())); +unwrapPhase(TContainerInOut &phase) { + const auto pi = std::numbers::pi_v; + auto prev = phase.front(); + std::transform(phase.begin() + 1, phase.end(), phase.begin() + 1, [&prev, pi](T ¤t) { + T diff = current - prev; + while (diff > pi) { + current -= 2 * pi; + diff = current - prev; + } + while (diff < -pi) { + current += 2 * pi; + diff = current - prev; + } + prev = current; + return current; + }); +} + +template TContainerOut = std::vector, + typename T = TContainerIn::value_type> + requires(std::is_same_v> || std::is_same_v>) +auto +computePhaseSpectrum(const TContainerIn &fftIn, TContainerOut &&phaseOut = {}, ConfigPhase config = {}) { + std::size_t phaseSize = config.computeHalfSpectrum ? fftIn.size() : (fftIn.size() / 2); + if constexpr (requires(std::size_t n) { phaseOut.resize(n); }) { + if (phaseOut.size() != phaseSize) { + phaseOut.resize(phaseSize); + } + } else { + static_assert(std::tuple_size_v == std::tuple_size_v, "Size mismatch for fixed-size container."); + } + std::transform(fftIn.begin(), std::next(fftIn.begin(), static_cast(phaseOut.size())), phaseOut.begin(), [](const auto &c) { return std::atan2(c.imag(), c.real()); }); + + if (config.unwrapPhase) { + unwrapPhase(phaseOut); + } + + if (config.outputInDeg) { + std::transform(phaseOut.begin(), phaseOut.end(), phaseOut.begin(), + [](const auto &phase) { return phase * static_cast(180.) * std::numbers::inv_pi_v; }); } - std::transform(fftOut.begin(), std::next(fftOut.begin(), static_cast(phaseSpectrum.size())), phaseSpectrum.begin(), - [](const auto &c) { return static_cast(std::atan2(c.imag(), c.real())); }); + + return phaseOut; +} + +template + requires(std::is_same_v> || std::is_same_v>) +auto +computePhaseSpectrum(const TContainerIn &fftIn, ConfigPhase config) { + return computePhaseSpectrum(fftIn, {}, config); } -} // namespace gr::algorithm +} // namespace gr::algorithm::fft #endif // GNURADIO_ALGORITHM_FFT_COMMON_HPP diff --git a/algorithm/include/gnuradio-4.0/algorithm/fourier/fft_types.hpp b/algorithm/include/gnuradio-4.0/algorithm/fourier/fft_types.hpp deleted file mode 100644 index b67367881..000000000 --- a/algorithm/include/gnuradio-4.0/algorithm/fourier/fft_types.hpp +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef GNURADIO_ALGORITHM_FFT_TYPES_HPP -#define GNURADIO_ALGORITHM_FFT_TYPES_HPP - -#include -#include - -namespace gr::algorithm { -template -concept ComplexType = std::is_same_v> || std::is_same_v>; - -template -using FFTInDataType = std::conditional_t, std::complex, U>; - -template -using FFTOutDataType = std::complex; - -template -struct FFTAlgoPrecision { - using type = T; -}; - -template -struct FFTAlgoPrecision { - using type = T::value_type; -}; - -} // namespace gr::algorithm -#endif // GNURADIO_ALGORITHM_FFT_TYPES_HPP diff --git a/algorithm/include/gnuradio-4.0/algorithm/fourier/fftw.hpp b/algorithm/include/gnuradio-4.0/algorithm/fourier/fftw.hpp index 66a9d058b..c4648e287 100644 --- a/algorithm/include/gnuradio-4.0/algorithm/fourier/fftw.hpp +++ b/algorithm/include/gnuradio-4.0/algorithm/fourier/fftw.hpp @@ -3,26 +3,22 @@ #include -#include "fft_types.hpp" #include "window.hpp" namespace gr::algorithm { -template -concept FFTwDoubleType = std::is_same_v> || std::is_same_v; - -template - requires(ComplexType || std::floating_point) +template, TInput, std::complex>> + requires((gr::meta::complex_like || std::floating_point) && (gr::meta::complex_like) ) struct FFTw { private: inline static std::mutex fftw_plan_mutex; public: // clang-format off - template - struct fftwImpl { + template + struct FFTwImpl { using PlanType = fftwf_plan; - using InAlgoDataType = std::conditional_t, fftwf_complex, float>; + using InAlgoDataType = std::conditional_t, fftwf_complex, float>; using OutAlgoDataType = fftwf_complex; using InUniquePtr = std::unique_ptr; using OutUniquePtr = std::unique_ptr; @@ -51,10 +47,11 @@ struct FFTw { static void forgetWisdom() {fftwf_forget_wisdom();} }; - template - struct fftwImpl { + template + requires (std::is_same_v> || std::is_same_v) + struct FFTwImpl { using PlanType = fftw_plan; - using InAlgoDataType = std::conditional_t, fftw_complex, double>; + using InAlgoDataType = std::conditional_t, fftw_complex, double>; using OutAlgoDataType = fftw_complex; using InUniquePtr = std::unique_ptr; using OutUniquePtr = std::unique_ptr; @@ -84,12 +81,13 @@ struct FFTw { }; // clang-format on - using OutDataType = std::conditional_t, T, std::complex>; - 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; + // Precision of the algorithm is defined by the output type `TOutput:value_type` + using AlgoDataType = std::conditional_t, TOutput, typename TOutput::value_type>; + 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" }; @@ -111,19 +109,16 @@ struct FFTw { ~FFTw() { clearFftw(); } - std::vector - computeFFT(const std::vector &in) { - if (fftSize != in.size()) { - fftSize = in.size(); - initAll(); + auto + compute(const std::ranges::input_range auto &in, std::ranges::output_range auto &&out) { + if constexpr (requires(std::size_t n) { out.resize(n); }) { + if (out.size() != in.size()) { + out.resize(in.size()); + } + } else { + static_assert(std::tuple_size_v == std::tuple_size_v, "Size mismatch for fixed-size container."); } - std::vector out(getOutputSize()); - computeFFT(in, out); - return out; - } - void - computeFFT(const std::vector &in, std::vector &out) { if (!std::has_single_bit(in.size())) { throw std::invalid_argument(fmt::format("Input data must have 2^N samples, input size: ", in.size())); } @@ -133,50 +128,69 @@ struct FFTw { initAll(); } - if (out.size() < getOutputSize()) { - throw std::out_of_range(fmt::format("Output vector size ({}) is not enough, at least {} needed. ", out.size(), getOutputSize())); + if (out.size() < fftSize) { + throw std::out_of_range(fmt::format("Output vector size ({}) is not enough, at least {} needed. ", out.size(), fftSize)); } - 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); + // precision is defined by output type, if needed cast input + if constexpr (!std::is_same_v) { + std::span inSpan(reinterpret_cast(fftwIn.get()), in.size()); + std::ranges::transform(in.begin(), in.end(), inSpan.begin(), [](const auto c) { return static_cast(c); }); + } else { + std::memcpy(fftwIn.get(), &(*in.begin()), sizeof(InAlgoDataType) * fftSize); + } + + FFTwImpl::execute(fftwPlan.get()); + + static_assert(sizeof(TOutput) == sizeof(OutAlgoDataType), "Sizes of TOutput type and OutAlgoDataType are not equal."); + std::memcpy(out.data(), fftwOut.get(), sizeof(TOutput) * getOutputSize()); - fftwImpl::execute(fftwPlan.get()); + // for the real input to complex a Hermitian output is produced by fftw, perform mirroring and conjugation fftw spectra to the second half + if (!gr::meta::complex_like) { + const auto halfIt = std::next(out.begin(), static_cast(fftSize / 2)); + std::ranges::transform(out.begin(), halfIt, halfIt, [](auto c) { return std::conj(c); }); + std::reverse(halfIt, out.end()); + } + + return out; + } - 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()); + auto + compute(const std::ranges::input_range auto &in) { + return compute(in, std::vector()); } [[nodiscard]] inline int importWisdom() const { // lock file while importing wisdom? - return fftwImpl::importWisdomFromFilename(wisdomPath); + return FFTwImpl::importWisdomFromFilename(wisdomPath); } [[nodiscard]] inline int exportWisdom() const { // lock file while exporting wisdom? - return fftwImpl::exportWisdomToFilename(wisdomPath); + return FFTwImpl::exportWisdomToFilename(wisdomPath); } [[nodiscard]] inline int importWisdomFromString(const std::string &wisdomString) const { - return fftwImpl::importWisdomFromString(wisdomString); + return FFTwImpl::importWisdomFromString(wisdomString); } [[nodiscard]] std::string exportWisdomToString() const { - return fftwImpl::exportWisdomToString(); + return FFTwImpl::exportWisdomToString(); } inline void forgetWisdom() const { - return fftwImpl::forgetWisdom(); + return FFTwImpl::forgetWisdom(); } private: [[nodiscard]] constexpr std::size_t getOutputSize() const { - if constexpr (ComplexType) { + if constexpr (gr::meta::complex_like) { return fftSize; } else { return 1 + fftSize / 2; @@ -186,14 +200,14 @@ struct FFTw { void initAll() { clearFftw(); - fftwIn = InUniquePtr(static_cast(fftwImpl::malloc(sizeof(InAlgoDataType) * fftSize))); - fftwOut = OutUniquePtr(static_cast(fftwImpl::malloc(sizeof(OutAlgoDataType) * getOutputSize()))); + 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 }; // what to do if error is returned std::ignore = importWisdom(); - fftwPlan = PlanUniquePtr(fftwImpl::plan(static_cast(fftSize), fftwIn.get(), fftwOut.get(), sign, flags)); + fftwPlan = PlanUniquePtr(FFTwImpl::plan(static_cast(fftSize), fftwIn.get(), fftwOut.get(), sign, flags)); std::ignore = exportWisdom(); } } diff --git a/algorithm/test/qa_algorithm_fourier.cpp b/algorithm/test/qa_algorithm_fourier.cpp index 3ac0e9b9d..c4bedbe33 100644 --- a/algorithm/test/qa_algorithm_fourier.cpp +++ b/algorithm/test/qa_algorithm_fourier.cpp @@ -1,6 +1,7 @@ #include +#include #include -#include +#include #include @@ -8,23 +9,186 @@ #include #include +#include +#include + +template +std::vector +generateSinSample(std::size_t N, double sample_rate, double frequency, double amplitude) { + std::vector signal(N); + for (std::size_t i = 0; i < N; i++) { + if constexpr (gr::meta::complex_like) { + signal[i] = { static_cast(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast(i) / sample_rate)), 0. }; + } else { + signal[i] = static_cast(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast(i) / sample_rate)); + } + } + return signal; +} template bool equalVectors(const T &v1, const U &v2, double tolerance = std::is_same_v ? 1.e-5 : 1e-4) { - if constexpr (gr::algorithm::ComplexType) { + if (v1.size() != v2.size()) { + return false; + } + if constexpr (gr::meta::complex_like) { return std::ranges::equal(v1, v2, [&tolerance](const auto &l, const auto &r) { - return std::abs(l.real() - r.real()) < static_cast(tolerance) && std::abs(l.imag() - r.imag()) < static_cast(tolerance); + return std::abs(l.real() - r.real()) < static_cast(tolerance) && std::abs(l.imag() - r.imag()) < static_cast(tolerance); }); } else { return std::ranges::equal(v1, v2, [&tolerance](const auto &l, const auto &r) { return std::abs(static_cast(l) - static_cast(r)) < tolerance; }); } } -const boost::ut::suite<"window functions"> windowTests = [] { +template +void +testFFTwTypes() { + using namespace boost::ut; + gr::algorithm::FFTw fftBlock; + expect(std::is_same_v, TExpInput>) << ""; + expect(std::is_same_v, TExpOutput>) << ""; + expect(std::is_same_v) << ""; +} + +template typename TAlgo> +struct TestTypes { + using InType = TInput; + using OutType = TOutput; + using AlgoType = TAlgo; +}; + +const boost::ut::suite<"FFT algorithms and window functions"> windowTests = [] { using namespace boost::ut; using namespace boost::ut::reflection; using gr::algorithm::window::create; + using gr::algorithm::FFT; + using gr::algorithm::FFTw; + + using ComplexTypesToTest = std::tuple< + // complex input, same in-out precision + TestTypes, std::complex, FFT>, TestTypes, std::complex, FFTw>, TestTypes, std::complex, FFT>, + TestTypes, std::complex, FFTw>, + // complex input, different in-out precision + TestTypes, std::complex, FFT>, TestTypes, std::complex, FFTw>, TestTypes, std::complex, FFT>, + TestTypes, std::complex, FFTw>>; + + using RealTypesToTest = std::tuple< + // real input, same in-out precision + TestTypes, FFT>, TestTypes, FFTw>, TestTypes, FFT>, TestTypes, FFTw>, + // real input, different in-out precision + TestTypes, FFT>, TestTypes, FFTw>, TestTypes, FFT>, TestTypes, FFTw>>; + + using AllTypesToTest = decltype(std::tuple_cat(std::declval(), std::declval())); + + "FFT algo sin tests"_test = []() { + typename T::AlgoType fftAlgo{}; + constexpr double tolerance{ 1.e-5 }; + struct TestParams { + std::uint32_t N{ 1024 }; // must be power of 2 + double sample_rate{ 128. }; // must be power of 2 (only for the unit test for easy comparison with true result) + double frequency{ 1. }; + double amplitude{ 1. }; + bool outputInDb{ false }; + }; + + std::vector testCases = { { 256, 128., 10., 5., false }, { 512, 4., 1., 1., false }, { 512, 32., 1., 0.1, false }, { 256, 128., 10., 5., false } }; + for (const auto &t : testCases) { + assert(std::has_single_bit(t.N)); + assert(std::has_single_bit(static_cast(t.sample_rate))); + + const auto signal{ generateSinSample(t.N, t.sample_rate, t.frequency, t.amplitude) }; + auto fftResult = fftAlgo.compute(signal); + auto magnitudeSpectrum = gr::algorithm::fft::computeMagnitudeSpectrum(fftResult); + auto phase = gr::algorithm::fft::computePhaseSpectrum(fftResult, { .outputInDeg = true, .unwrapPhase = true }); + const auto peakIndex{ static_cast( + std::distance(magnitudeSpectrum.begin(), + std::max_element(magnitudeSpectrum.begin(), std::next(magnitudeSpectrum.begin(), static_cast(t.N / 2u))))) }; // only positive frequencies from FFT + const auto peakAmplitude = magnitudeSpectrum[peakIndex]; + const auto peakFrequency{ static_cast(peakIndex) * t.sample_rate / static_cast(t.N) }; + + const auto expectedAmplitude = t.outputInDb ? 20. * log10(std::abs(t.amplitude)) : t.amplitude; + 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()); + } + } | AllTypesToTest{}; + + "FFT algo pattern tests"_test = []() { + using InType = T::InType; + typename T::AlgoType fftAlgo{}; + constexpr double tolerance{ 1.e-5 }; + constexpr std::uint32_t N{ 16 }; + static_assert(N == 16, "expected values are calculated for N == 16"); + + std::vector signal(N); + std::size_t expectedPeakIndex{ 0 }; + InType expectedFft0{ 0., 0. }; + double expectedPeakAmplitude{ 0. }; + for (std::size_t iT = 0; iT < 5; iT++) { + if (iT == 0) { + std::ranges::fill(signal.begin(), signal.end(), InType(0., 0.)); + expectedFft0 = { 0., 0. }; + expectedPeakAmplitude = 0.; + } else if (iT == 1) { + std::ranges::fill(signal.begin(), signal.end(), InType(1., 0.)); + expectedFft0 = { 16., 0. }; + expectedPeakAmplitude = 2.; + } else if (iT == 2) { + std::ranges::fill(signal.begin(), signal.end(), InType(1., 1.)); + expectedFft0 = { 16., 16. }; + expectedPeakAmplitude = std::sqrt(8.); + } else if (iT == 3) { + std::iota(signal.begin(), signal.end(), 1); + expectedFft0 = { 136., 0. }; + expectedPeakAmplitude = 17.; + } else if (iT == 4) { + int i = 0; + std::ranges::generate(signal.begin(), signal.end(), [&i] { return InType(static_cast(i++ % 2), 0.); }); + expectedFft0 = { 8., 0. }; + expectedPeakAmplitude = 1.; + } + + auto fftResult = fftAlgo.compute(signal); + auto magnitudeSpectrum = gr::algorithm::fft::computeMagnitudeSpectrum(fftResult); + + const auto peakIndex{ static_cast(std::distance(magnitudeSpectrum.begin(), std::ranges::max_element(magnitudeSpectrum))) }; + const auto peakAmplitude{ magnitudeSpectrum[peakIndex] }; + + 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(fftResult[0].real()), static_cast(expectedFft0.real()), tolerance)) << fmt::format("<{}> equal fft[0].real()", type_name()); + expect(approx(static_cast(fftResult[0].imag()), static_cast(expectedFft0.imag()), tolerance)) << fmt::format("<{}> equal fft[0].imag()", type_name()); + } + } | ComplexTypesToTest{}; + + "Unwrap Phase tests"_test = [] { + std::vector phase = { 0.2, -1., 2.5, -3.1, 0.9, -0.5, 1.2, 0.8, 1.5, -1.2, -2.7, 0.9, -0.8, -1.4, 0.6, 1.1, -1.9, 0.4, 1.3, -0.7 }; + // Output generated with python numpy.unwrap(phase) + std::vector expOut = { 0.2, -1., -3.78318531, -3.1, -5.38318531, -6.78318531, -5.08318531, -5.48318531, -4.78318531, -7.48318531, + -8.98318531, -11.66637061, -13.36637061, -13.96637061, -11.96637061, -11.46637061, -14.46637061, -12.16637061, -11.26637061, -13.26637061 }; + gr::algorithm::fft::unwrapPhase(phase); + expect(equalVectors(phase, expOut)) << "unwrapped phases are equal"; + }; + + "FFTw types tests"_test = [] { + testFFTwTypes, std::complex, fftwf_complex, fftwf_complex, fftwf_plan>(); + testFFTwTypes, std::complex, fftw_complex, fftw_complex, fftw_plan>(); + testFFTwTypes, float, fftwf_complex, fftwf_plan>(); + testFFTwTypes, double, fftw_complex, fftw_plan>(); + }; + + "FFTW wisdom import/export tests"_test = []() { + gr::algorithm::FFTw> fftw1{}; + + std::string wisdomString1 = fftw1.exportWisdomToString(); + fftw1.forgetWisdom(); + int importOk = fftw1.importWisdomFromString(wisdomString1); + expect(eq(importOk, 1)) << "Wisdom import from string."; + 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."; + }; "window pre-computed array tests"_test = []() { // this tests regression w.r.t. changed implementations // Expected value for size 8 diff --git a/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp b/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp index ba8082491..7f67c1613 100644 --- a/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp +++ b/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp @@ -9,7 +9,6 @@ #include #include -#include #include #include @@ -17,8 +16,18 @@ namespace gr::blocks::fft { using namespace gr; -template, typename FourierAlgorithm = gr::algorithm::FFT>> - requires(gr::algorithm::ComplexType || std::floating_point || std::is_same_v> || std::is_same_v>) +template +struct OutputDataSet { + using type = DataSet; +}; + +template +struct OutputDataSet { + using type = DataSet; +}; + +template::type, template typename FourierAlgorithm = gr::algorithm::FFT> + requires((gr::meta::complex_like || std::floating_point) && (std::is_same_v> || std::is_same_v>) ) struct FFT : public Block, ResamplingRatio<1LU, 1024LU>, Doc> { using value_type = U::value_type; - using InDataType = gr::algorithm::FFTInDataType; - using OutDataType = gr::algorithm::FFTOutDataType; + using InDataType = std::conditional_t, std::complex, value_type>; + using OutDataType = std::complex; - PortIn in; - PortOut out; + PortIn in; + PortOut out; - FourierAlgorithm _fftImpl; - gr::algorithm::window::Type _windowType = gr::algorithm::window::Type::Hann; - std::vector _window = gr::algorithm::window::create(_windowType, 1024U); + FourierAlgorithm> _fftImpl; + gr::algorithm::window::Type _windowType = gr::algorithm::window::Type::Hann; + std::vector _window = gr::algorithm::window::create(_windowType, 1024U); // settings - const std::string algorithm = gr::meta::type_name(); + const std::string algorithm = gr::meta::type_name(); Annotated> fftSize{ 1024U }; Annotated> window = std::string(gr::algorithm::window::to_string(_windowType)); Annotated> outputInDb{ false }; + Annotated> outputInDeg{ false }; + Annotated> unwrapPhase{ false }; Annotated, Unit<"Hz">> sample_rate = 1.f; Annotated signal_name = "unknown signal"; Annotated> signal_unit = "a.u."; @@ -101,10 +112,11 @@ On the choice of window (mathematically aka. apodisation) functions Annotated> signal_max = std::numeric_limits::max(); // semi-private caching vectors (need to be public for unit-test) -> TODO: move to FFT implementations, casting from T -> U::value_type should be done there - std::vector _inData = std::vector(fftSize, 0); - std::vector _outData = std::vector(gr::algorithm::ComplexType ? fftSize.value : (1U + fftSize.value / 2U), 0); - std::vector _magnitudeSpectrum = std::vector(_outData.size(), 0); - std::vector _phaseSpectrum = std::vector(_outData.size(), 0); + std::vector _inData = std::vector(fftSize, 0); + std::vector _outData = std::vector(gr::meta::complex_like ? fftSize.value : (1U + fftSize.value / 2U), 0); + std::vector _magnitudeSpectrum = std::vector(_outData.size(), 0); + std::vector _phaseSpectrum = std::vector(_outData.size(), 0); + constexpr static bool computeHalfSpectrum = gr::meta::complex_like; void settingsChanged(const property_map & /*old_settings*/, const property_map &newSettings) noexcept { @@ -124,7 +136,6 @@ On the choice of window (mathematically aka. apodisation) functions // N.B. this should become part of the Fourier transform implementation _inData.resize(fftSize, 0); - constexpr bool computeHalfSpectrum = gr::algorithm::ComplexType; _outData.resize(computeHalfSpectrum ? newSize : (1U + newSize / 2), 0); _magnitudeSpectrum.resize(computeHalfSpectrum ? newSize : (newSize / 2), 0); _phaseSpectrum.resize(computeHalfSpectrum ? newSize : (newSize / 2), 0); @@ -140,7 +151,7 @@ On the choice of window (mathematically aka. apodisation) functions // apply window function for (std::size_t i = 0U; i < fftSize; i++) { - if constexpr (gr::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { _inData[i].real(_inData[i].real() * _window[i]); _inData[i].imag(_inData[i].imag() * _window[i]); } else { @@ -148,15 +159,13 @@ On the choice of window (mathematically aka. apodisation) functions } } - _outData = _fftImpl.computeFFT(_inData); + _outData = _fftImpl.compute(_inData); + _magnitudeSpectrum = gr::algorithm::fft::computeMagnitudeSpectrum(_outData, _magnitudeSpectrum, + algorithm::fft::ConfigMagnitude{ .computeHalfSpectrum = computeHalfSpectrum, .outputInDb = outputInDb }); + _phaseSpectrum = gr::algorithm::fft::computePhaseSpectrum(_outData, _phaseSpectrum, + algorithm::fft::ConfigPhase{ .computeHalfSpectrum = computeHalfSpectrum, .outputInDeg = outputInDeg, .unwrapPhase = unwrapPhase }); - gr::algorithm::computeMagnitudeSpectrum(_outData, _magnitudeSpectrum, fftSize, outputInDb); - gr::algorithm::computePhaseSpectrum(_outData, _phaseSpectrum); - if constexpr (std::is_same_v> || std::is_same_v>) { - output[0] = createDataset(); - } else { - static_assert(!std::is_same_v> && "FFT output type not (yet) implemented"); - } + output[0] = createDataset(); return work::Status::OK; } @@ -178,12 +187,12 @@ On the choice of window (mathematically aka. apodisation) functions ds.signal_values.resize(dim * N); auto const freqWidth = static_cast(sample_rate) / static_cast(fftSize); - if constexpr (gr::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { auto const freqOffset = static_cast(N / 2) * freqWidth; std::ranges::transform(std::views::iota(0UL, N), std::ranges::begin(ds.signal_values), [freqWidth, freqOffset](const auto i) { return static_cast(i) * freqWidth - freqOffset; }); } else { - std::ranges::transform(std::views::iota(0UL, N), std::ranges::begin(ds.signal_values), [freqWidth](const auto i) { return static_cast(i) * freqWidth; }); + std::ranges::transform(std::views::iota(0UL, N), std::ranges::begin(ds.signal_values), [freqWidth](const auto i) { return static_cast(i) * freqWidth; }); } std::ranges::transform(_outData.begin(), _outData.end(), std::next(ds.signal_values.begin(), static_cast(N)), [](const auto &c) { return c.real(); }); std::ranges::transform(_outData.begin(), _outData.end(), std::next(ds.signal_values.begin(), static_cast(2U * N)), [](const auto &c) { return c.imag(); }); @@ -205,6 +214,8 @@ On the choice of window (mathematically aka. apodisation) functions { "fft_size", fftSize }, { "window", window }, { "output_in_db", outputInDb }, + { "output_in_deg", outputInDeg }, + { "unwrap_phase", unwrapPhase }, { "numerator", this->numerator }, { "denominator", this->denominator }, { "stride", this->stride } } }; @@ -215,7 +226,7 @@ On the choice of window (mathematically aka. apodisation) functions } // namespace gr::blocks::fft -ENABLE_REFLECTION_FOR_TEMPLATE_FULL((typename T, typename U, typename FourierAlgoImpl), (gr::blocks::fft::FFT), // - in, out, algorithm, fftSize, window, outputInDb, sample_rate, signal_name, signal_unit, signal_min, signal_max); +ENABLE_REFLECTION_FOR_TEMPLATE_FULL((typename T, typename U, template typename FourierAlgoImpl), (gr::blocks::fft::FFT), // + in, out, algorithm, fftSize, window, outputInDb, outputInDeg, unwrapPhase, sample_rate, signal_name, signal_unit, signal_min, signal_max); #endif // GNURADIO_FFT_HPP diff --git a/blocks/fourier/test/qa_fourier.cpp b/blocks/fourier/test/qa_fourier.cpp index cc8fe7a6e..7d958eba4 100644 --- a/blocks/fourier/test/qa_fourier.cpp +++ b/blocks/fourier/test/qa_fourier.cpp @@ -44,7 +44,7 @@ std::vector generateSinSample(std::size_t N, double sample_rate, double frequency, double amplitude) { std::vector signal(N); for (std::size_t i = 0; i < N; i++) { - if constexpr (gr::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { signal[i] = { static_cast(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast(i) / sample_rate)), 0. }; } else { signal[i] = static_cast(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast(i) / sample_rate)); @@ -59,7 +59,7 @@ equalVectors(const std::vector &v1, const std::vector &v2, double toleranc if (v1.size() != v2.size()) { return false; } - if constexpr (gr::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { return std::equal(v1.begin(), v1.end(), v2.begin(), [&tolerance](const auto &l, const auto &r) { return std::abs(l.real() - r.real()) < static_cast(tolerance) && std::abs(l.imag() - r.imag()) < static_cast(tolerance); }); @@ -68,19 +68,9 @@ equalVectors(const std::vector &v1, const std::vector &v2, double toleranc } } -template +template void -testFFTwTypes() { - using namespace boost::ut; - gr::algorithm::FFTw fftBlock; - 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, A> &fftBlock, const gr::DataSet &ds1, float sample_rate) { +equalDataset(const gr::blocks::fft::FFT> &fftBlock, const gr::DataSet &ds1, float sample_rate) { using namespace boost::ut; using namespace boost::ut::reflection; @@ -118,17 +108,10 @@ equalDataset(const gr::blocks::fft::FFT, A> &fftBlock, const g } } -template -using FFTwAlgo = gr::algorithm::FFTw; - -template -using FFTAlgo = gr::algorithm::FFT; - -template typename A> -struct TypePair { - using InType = T; - using OutType = U; - using AlgoType = A>; +template +struct TestTypes { + using InType = TInput; + using OutType = TOutput; }; const boost::ut::suite<"Fourier Transforms"> fftTests = [] { @@ -136,116 +119,26 @@ const boost::ut::suite<"Fourier Transforms"> fftTests = [] { using namespace gr::blocks::fft; using namespace boost::ut::reflection; - std::tuple, DataSet, FFTwAlgo>, TypePair, DataSet, FFTAlgo>, TypePair, DataSet, FFTwAlgo>, - TypePair, DataSet, FFTAlgo>> - complexTypesWithAlgoToTest{}; - - std::tuple, DataSet, FFTwAlgo>, TypePair, DataSet, FFTAlgo>, TypePair, DataSet, FFTwAlgo>, - TypePair, DataSet, FFTAlgo>, TypePair, FFTwAlgo>, TypePair, FFTAlgo>, TypePair, FFTwAlgo>, - TypePair, FFTAlgo>> - typesWithAlgoToTest{}; - - "FFT sin tests"_test = []() { - using InType = T::InType; - using OutType = T::OutType; - using AlgoType = T::AlgoType; - FFT fftBlock{}; - constexpr double tolerance{ 1.e-5 }; - struct TestParams { - std::uint32_t N{ 1024 }; // must be power of 2 - double sample_rate{ 128. }; // must be power of 2 (only for the unit test for easy comparison with true result) - double frequency{ 1. }; - double amplitude{ 1. }; - bool outputInDb{ false }; - }; - - std::vector testCases = { { 256, 128., 10., 5., false }, { 512, 4., 1., 1., false }, { 512, 32., 1., 0.1, false }, { 256, 128., 10., 5., true } }; - for (const auto &t : testCases) { - assert(std::has_single_bit(t.N)); - assert(std::has_single_bit(static_cast(t.sample_rate))); - - std::ignore = fftBlock.settings().set({ { "fftSize", t.N } }); - std::ignore = fftBlock.settings().set({ { "outputInDb", t.outputInDb } }); - std::ignore = fftBlock.settings().set({ { "window", "None" } }); - std::ignore = fftBlock.settings().applyStagedParameters(); - const auto signal{ generateSinSample(t.N, t.sample_rate, t.frequency, t.amplitude) }; - std::vector resultingDataSets(1); - expect(gr::work::Status::OK == fftBlock.processBulk(signal, resultingDataSets)); - - const auto peakIndex{ - static_cast(std::distance(fftBlock._magnitudeSpectrum.begin(), - std::max_element(fftBlock._magnitudeSpectrum.begin(), std::next(fftBlock._magnitudeSpectrum.begin(), static_cast(t.N / 2u))))) - }; // only positive frequencies from FFT - const auto peakAmplitude = fftBlock._magnitudeSpectrum[peakIndex]; - const auto peakFrequency{ static_cast(peakIndex) * t.sample_rate / static_cast(t.N) }; - - const auto expectedAmplitude = t.outputInDb ? 20. * log10(std::abs(t.amplitude)) : t.amplitude; - 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()); - } - } | typesWithAlgoToTest; - - "FFT pattern tests"_test = []() { - using InType = T::InType; - using OutType = T::OutType; - using AlgoType = T::AlgoType; - constexpr double tolerance{ 1.e-5 }; - constexpr std::uint32_t N{ 16 }; - FFT fftBlock({ { "fftSize", N }, { "window", "None" } }); - std::ignore = fftBlock.settings().applyStagedParameters(); - - std::vector signal(N); - - static_assert(N == 16, "expected values are calculated for N == 16"); - std::size_t expectedPeakIndex{ 0 }; - InType expectedFft0{ 0., 0. }; - double expectedPeakAmplitude{ 0. }; - for (std::size_t iT = 0; iT < 5; iT++) { - if (iT == 0) { - std::ranges::fill(signal.begin(), signal.end(), InType(0., 0.)); - expectedFft0 = { 0., 0. }; - expectedPeakAmplitude = 0.; - } else if (iT == 1) { - std::ranges::fill(signal.begin(), signal.end(), InType(1., 0.)); - expectedFft0 = { 16., 0. }; - expectedPeakAmplitude = 2.; - } else if (iT == 2) { - std::ranges::fill(signal.begin(), signal.end(), InType(1., 1.)); - expectedFft0 = { 16., 16. }; - expectedPeakAmplitude = std::sqrt(8.); - } else if (iT == 3) { - std::iota(signal.begin(), signal.end(), 1); - expectedFft0 = { 136., 0. }; - expectedPeakAmplitude = 17.; - } else if (iT == 4) { - int i = 0; - std::ranges::generate(signal.begin(), signal.end(), [&i] { return InType(static_cast(i++ % 2), 0.); }); - expectedFft0 = { 8., 0. }; - expectedPeakAmplitude = 1.; - } - std::vector resultingDataSets(1); - expect(gr::work::Status::OK == fftBlock.processBulk(signal, resultingDataSets)); - - const auto peakIndex{ static_cast(std::distance(fftBlock._magnitudeSpectrum.begin(), std::ranges::max_element(fftBlock._magnitudeSpectrum))) }; - const auto peakAmplitude{ fftBlock._magnitudeSpectrum[peakIndex] }; - - 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(fftBlock._outData[0].real()), static_cast(expectedFft0.real()), tolerance)) << fmt::format("<{}> equal fft[0].real()", type_name()); - expect(approx(static_cast(fftBlock._outData[0].imag()), static_cast(expectedFft0.imag()), tolerance)) << fmt::format("<{}> equal fft[0].imag()", type_name()); - } - } | complexTypesWithAlgoToTest; + using AllTypesToTest = std::tuple< + // complex input, same in-out precision + TestTypes, DataSet>, TestTypes, DataSet>, + // complex input, different in-out precision + TestTypes, DataSet>, TestTypes, DataSet>, + // real input, same in-out precision + TestTypes>, TestTypes>, + // real input, different in-out precision + TestTypes>, TestTypes>>; "FFT processBulk tests"_test = []() { - using InType = T::InType; - using OutType = T::OutType; - using AlgoType = T::AlgoType; + using InType = T::InType; + using OutType = T::OutType; - constexpr std::uint32_t N{ 16 }; - constexpr float sample_rate{ 1.f }; - FFT fftBlock({ { "fftSize", N }, { "sample_rate", sample_rate } }); + constexpr std::uint32_t N{ 16 }; + constexpr float sample_rate{ 1.f }; + FFT fftBlock({ { "fftSize", N }, { "sample_rate", sample_rate } }); std::ignore = fftBlock.settings().applyStagedParameters(); - expect(eq(fftBlock.algorithm, gr::meta::type_name())); + + expect(eq(fftBlock.algorithm, gr::meta::type_name>>())); std::vector signal(N); std::iota(signal.begin(), signal.end(), 1); @@ -254,68 +147,47 @@ const boost::ut::suite<"Fourier Transforms"> fftTests = [] { expect(gr::work::Status::OK == fftBlock.processBulk(signal, outSpan)); equalDataset(fftBlock, v[0], sample_rate); - } | typesWithAlgoToTest; + } | AllTypesToTest{}; - "FFT types tests"_test = [] { + "FFT block types tests"_test = [] { expect(std::is_same_v>::value_type, float>) << "output type must be float"; - expect(std::is_same_v>::value_type, float>) << "output type must be float"; + expect(std::is_same_v>::value_type, double>) << "output type must be float"; expect(std::is_same_v::value_type, float>) << "output type must be float"; - expect(std::is_same_v::value_type, float>) << "output type must be float"; - expect(std::is_same_v::value_type, float>) << "output type must be float"; + expect(std::is_same_v::value_type, double>) << "output type must be float"; expect(std::is_same_v, gr::DataSet>::value_type, double>) << "output type must be double"; expect(std::is_same_v>::value_type, double>) << "output type must be double"; }; - "FFT fftw types tests"_test = [] { - testFFTwTypes, fftwf_complex, fftwf_complex, fftwf_plan>(); - testFFTwTypes, fftw_complex, fftw_complex, fftw_plan>(); - testFFTwTypes(); - testFFTwTypes(); - }; - "FFT flow graph example"_test = [] { // This test checks how fftw works if one creates and destroys several fft blocks in different graph flows using namespace boost::ut; using Scheduler = gr::scheduler::Simple<>; auto threadPool = std::make_shared("custom pool", gr::thread_pool::CPU_BOUND, 2, 2); gr::Graph flow1; - auto &source1 = flow1.emplaceBlock>(); - auto &fftBlock = flow1.emplaceBlock>({ { "fftSize", static_cast(16) } }); + auto &source1 = flow1.emplaceBlock>(); + auto &fftBlock = flow1.emplaceBlock>({ { "fftSize", static_cast(16) } }); expect(eq(gr::ConnectionResult::SUCCESS, flow1.connect<"out">(source1).to<"in">(fftBlock))); auto sched1 = Scheduler(std::move(flow1), threadPool); // run 2 times to check potential memory problems for (int i = 0; i < 2; i++) { gr::Graph flow2; - auto &source2 = flow2.emplaceBlock>(); - auto &fft2 = flow2.emplaceBlock>({ { "fftSize", static_cast(16) } }); + auto &source2 = flow2.emplaceBlock>(); + auto &fft2 = flow2.emplaceBlock>({ { "fftSize", static_cast(16) } }); expect(eq(gr::ConnectionResult::SUCCESS, flow2.connect<"out">(source2).to<"in">(fft2))); auto sched2 = Scheduler(std::move(flow2), threadPool); sched2.runAndWait(); + expect(approx(source2.count, source2.nSamples, 1e-4)); } sched1.runAndWait(); expect(approx(source1.count, source1.nSamples, 1e-4)); }; - "FFTW wisdom import/export tests"_test = []() { - gr::algorithm::FFTw fftw1{}; - - std::string wisdomString1 = fftw1.exportWisdomToString(); - fftw1.forgetWisdom(); - int importOk = fftw1.importWisdomFromString(wisdomString1); - expect(eq(importOk, 1)) << "Wisdom import from string."; - 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."; - }; - "window function tests"_test = []() { - using InType = T::InType; - using OutType = T::OutType; - using AlgoType = T::AlgoType; + using InType = T::InType; + using OutType = T::OutType; - FFT fftBlock{}; + FFT fftBlock{}; using value_type = OutType::value_type; constexpr value_type tolerance{ value_type(0.00001) }; @@ -327,7 +199,7 @@ const boost::ut::suite<"Fourier Transforms"> fftTests = [] { std::ignore = fftBlock.settings().applyStagedParameters(); std::vector signal(N); - if constexpr (gr::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { typename InType::value_type i = 0.; std::ranges::generate(signal.begin(), signal.end(), [&i] { i = i + static_cast(1.); @@ -345,7 +217,7 @@ const boost::ut::suite<"Fourier Transforms"> fftTests = [] { std::vector windowFunc = gr::algorithm::window::create(window, N); for (std::size_t i = 0; i < N; i++) { - if constexpr (gr::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { const auto expValue = static_cast(signal[i].real()) * windowFunc[i]; expect(approx(fftBlock._inData[i].real(), expValue, tolerance)) << fmt::format("<{}> equal fftwIn complex.real", type_name()); expect(approx(fftBlock._inData[i].imag(), expValue, tolerance)) << fmt::format("<{}> equal fftwIn complex.imag", type_name()); @@ -355,7 +227,7 @@ const boost::ut::suite<"Fourier Transforms"> fftTests = [] { } } } - } | typesWithAlgoToTest; + } | AllTypesToTest{}; }; int diff --git a/core/benchmarks/bm_fft.cpp b/core/benchmarks/bm_fft.cpp index 3ccba915e..5fd10c5b5 100644 --- a/core/benchmarks/bm_fft.cpp +++ b/core/benchmarks/bm_fft.cpp @@ -17,9 +17,9 @@ * reference: https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#Pseudocode */ template - requires gr::algorithm::ComplexType + requires gr::meta::complex_like void -computeFFFTCooleyTukey(std::vector &signal) { +computeFFTCooleyTukey(std::vector &signal) { const std::size_t N{ signal.size() }; if (N == 1) return; @@ -31,8 +31,8 @@ computeFFFTCooleyTukey(std::vector &signal) { odd[i] = signal[2 * i + 1]; } - computeFFFTCooleyTukey(even); - computeFFFTCooleyTukey(odd); + computeFFTCooleyTukey(even); + computeFFTCooleyTukey(odd); 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++) { @@ -47,7 +47,7 @@ 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::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { 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)); @@ -56,6 +56,16 @@ generateSinSample(std::size_t N, double sampleRate, double frequency, double amp return signal; } +template +struct FFTAlgoPrecision { + using type = T; +}; + +template +struct FFTAlgoPrecision { + using type = T::value_type; +}; + template void testFFT() { @@ -78,7 +88,7 @@ testFFT() { std::vector signal = generateSinSample(N, sampleRate, frequency, amplitude); { - gr::blocks::fft::FFT, FFTw>> fft1({ { "fftSize", N } }); + gr::blocks::fft::FFT, FFTw> fft1({ { "fftSize", N } }); std::ignore = fft1.settings().applyStagedParameters(); std::vector> resultingDataSets(1); @@ -87,7 +97,7 @@ testFFT() { }; } { - gr::blocks::fft::FFT, FFT>> fft1({ { "fftSize", N } }); + gr::blocks::fft::FFT, FFT> fft1({ { "fftSize", N } }); std::ignore = fft1.settings().applyStagedParameters(); std::vector> resultingDataSets(1); @@ -96,10 +106,10 @@ testFFT() { }; } - if constexpr (gr::algorithm::ComplexType) { + if constexpr (gr::meta::complex_like) { ::benchmark::benchmark(fmt::format("{} - fftCT", type_name())) = [&signal] { auto signalCopy = signal; - computeFFFTCooleyTukey(signalCopy); + computeFFTCooleyTukey(signalCopy); }; } diff --git a/meta/include/gnuradio-4.0/meta/utils.hpp b/meta/include/gnuradio-4.0/meta/utils.hpp index 8639b2a7b..2c1cd0cdb 100644 --- a/meta/include/gnuradio-4.0/meta/utils.hpp +++ b/meta/include/gnuradio-4.0/meta/utils.hpp @@ -1,6 +1,7 @@ #ifndef GNURADIO_GRAPH_UTILS_HPP #define GNURADIO_GRAPH_UTILS_HPP +#include #include #include #include @@ -251,6 +252,9 @@ concept vectorizable_v = std::constructible_from>; template using vectorizable = std::integral_constant>; +template +concept complex_like = std::is_same_v> || std::is_same_v>; + /** * Determines the SIMD width of the given structure. This can either be a stdx::simd object or a * tuple-like of stdx::simd (recursively). The latter requires that the SIMD width is homogeneous.