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.
fair-acc · Sep 27, 2023 · 32fdfcd · 32fdfcd
1 parent ed71b4c
commit 32fdfcd
Show file tree

Hide file tree

Showing 10 changed files with 754 additions and 478 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -65,11 +65,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
@@ -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 <fmt/format.h>
 #include <numbers>
 
-/// 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<typename T>
-    requires gr::blocks::fft::ComplexType<T>
-void
-computeFftV1(std::vector<T> &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::size_t>(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<typename T::value_type>(-2. * std::numbers::pi) / static_cast<typename T::value_type>(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<typename T>
-    requires gr::blocks::fft::ComplexType<T>
+    requires gr::algorithm::ComplexType<T>
 void
-computeFftV2(std::vector<T> &signal) {
+computeFFFTCooleyTukey(std::vector<T> &signal) {
     const std::size_t N{ signal.size() };
 
     if (N == 1) return;
@@ -59,35 +25,23 @@ computeFftV2(std::vector<T> &signal) {
         odd[i]  = signal[2 * i + 1];
     }
 
-    computeFftV2(even);
-    computeFftV2(odd);
+    computeFFFTCooleyTukey(even);
+    computeFFFTCooleyTukey(odd);
 
-    const typename T::value_type wn{ static_cast<typename T::value_type>(2. * std::numbers::pi) / static_cast<typename T::value_type>(N) };
+    const typename T::value_type wn{ static_cast<typename T::value_type>(2. * std::numbers::pi_v<double>) / static_cast<typename T::value_type>(N) };
     for (std::size_t i = 0; i < N / 2; i++) {
         const T wkn(std::cos(wn * static_cast<typename T::value_type>(i)), std::sin(wn * static_cast<typename T::value_type>(i)));
         signal[i]         = even[i] + wkn * odd[i];
         signal[i + N / 2] = even[i] - wkn * odd[i];
     }
 }
 
-template<typename T>
-    requires gr::blocks::fft::ComplexType<T>
-std::vector<typename T::value_type>
-computeMagnitudeSpectrum(std::vector<T> &fftSignal) {
-    const std::size_t                   N{ fftSignal.size() };
-    std::vector<typename T::value_type> 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<typename T::value_type>(2.) / static_cast<typename T::value_type>(N);
-    }
-    return magnitudeSpectrum;
-}
-
 template<typename T>
 std::vector<T>
 generateSinSample(std::size_t N, double sampleRate, double frequency, double amplitude) {
     std::vector<T> signal(N);
     for (std::size_t i = 0; i < N; i++) {
-        if constexpr (gr::blocks::fft::ComplexType<T>) {
+        if constexpr (gr::algorithm::ComplexType<T>) {
             signal[i] = { static_cast<typename T::value_type>(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast<double>(i) / sampleRate)), 0. };
         } else {
             signal[i] = static_cast<T>(amplitude * std::sin(2. * std::numbers::pi * frequency * static_cast<double>(i) / sampleRate));
@@ -98,7 +52,7 @@ generateSinSample(std::size_t N, double sampleRate, double frequency, double amp
 
 template<typename T>
 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<T> signal  = generateSinSample<T>(N, sampleRate, frequency, amplitude);
-    std::string    nameOpt = withMagSpectrum ? "fft+mag" : "fft";
+    std::vector<T> signal = generateSinSample<T>(N, sampleRate, frequency, amplitude);
 
     {
-        gr::blocks::fft::fft<T> fft1{};
+        gr::blocks::fft::FFT<T, gr::algorithm::FFTw<T>> 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<nRepetitions>(fmt::format("{} - {} fftw", nameOpt, type_name<T>())) = [&fft1, &withMagSpectrum] {
+        ::benchmark::benchmark<nRepetitions>(fmt::format("{} - fftw", type_name<T>())) = [&fft1] {
+            fft1.prepareInput();
+            fft1.computeFFT();
+        };
+    }
+    {
+        gr::blocks::fft::FFT<T, gr::algorithm::FFT<T>> 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<nRepetitions>(fmt::format("{} - fft", type_name<T>())) = [&fft1] {
             fft1.prepareInput();
-            fft1.computeFft();
-            if (withMagSpectrum) fft1.computeMagnitudeSpectrum();
+            fft1.computeFFT();
         };
     }
 
     if constexpr (gr::blocks::fft::ComplexType<T>) {
-        ::benchmark::benchmark<nRepetitions>(fmt::format("{} - {} fft_v1", nameOpt, type_name<T>())) = [&signal, &withMagSpectrum] {
+        ::benchmark::benchmark<nRepetitions>(fmt::format("{} - fftCT", type_name<T>())) = [&signal] {
             auto signalCopy = signal;
-            computeFftV1<T>(signalCopy);
-            if (withMagSpectrum) [[maybe_unused]]
-                auto magnitudeSpectrum = computeMagnitudeSpectrum<T>(signalCopy);
+            computeFFFTCooleyTukey<T>(signalCopy);
         };
+    }
 
-        ::benchmark::benchmark<nRepetitions>(fmt::format("{} - {} fft_v2", nameOpt, type_name<T>())) = [&signal, &withMagSpectrum] {
-            auto signalCopy = signal;
-            computeFftV2<T>(signalCopy);
-            if (withMagSpectrum) [[maybe_unused]]
-                auto magnitudeSpectrum = computeMagnitudeSpectrum<T>(signalCopy);
+    {
+        gr::blocks::fft::FFT<T, gr::algorithm::FFTw<T>> 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<nRepetitions>(fmt::format("{} - magnitude", type_name<T>())) = [&fft1] { fft1.computeMagnitudeSpectrum(); };
+        ::benchmark::benchmark<nRepetitions>(fmt::format("{} - phase", type_name<T>()))     = [&fft1] { fft1.computePhaseSpectrum(); };
+
+        double sum{ 0. };
+        fft1.computeMagnitudeSpectrum();
+        fft1.computePhaseSpectrum();
+        ::benchmark::benchmark<nRepetitions>(fmt::format("{} - dataset", type_name<T>())) = [&fft1, &sum] {
+            const auto ds1 = fft1.createDataset();
+            sum += static_cast<double>(ds1.signal_values[0]);
         };
     }
+
+    ::benchmark::results::add_separator();
 }
 
 inline const boost::ut::suite _fft_bm_tests = [] {
     std::tuple<std::complex<float>, std::complex<double>> complexTypesToTest{};
     std::tuple<float, double>                             realTypesToTest{};
 
-    std::apply([]<class... TArgs>(TArgs... /*args*/) { (testFft<TArgs>(false), ...); }, complexTypesToTest);
-    benchmark::results::add_separator();
-    std::apply([]<class... TArgs>(TArgs... /*args*/) { (testFft<TArgs>(true), ...); }, complexTypesToTest);
-    benchmark::results::add_separator();
-    std::apply([]<class... TArgs>(TArgs... /*args*/) { (testFft<TArgs>(false), ...); }, realTypesToTest);
-    benchmark::results::add_separator();
-    std::apply([]<class... TArgs>(TArgs... /*args*/) { (testFft<TArgs>(true), ...); }, realTypesToTest);
+    std::apply([]<class... TArgs>(TArgs... /*args*/) { (testFFT<TArgs>(), ...); }, complexTypesToTest);
+    // benchmark::results::add_separator();
+    std::apply([]<class... TArgs>(TArgs... /*args*/) { (testFFT<TArgs>(), ...); }, realTypesToTest);
 };
 
 int

diff --git a/cmake/CompilerWarnings.cmake b/cmake/CompilerWarnings.cmake
@@ -3,83 +3,83 @@
 # https://github.com/lefticus/cppbestpractices/blob/master/02-Use_the_Tools_Available.md
 
 function(set_project_warnings project_name)
-  option(WARNINGS_AS_ERRORS "Treat compiler warnings as errors" FALSE) #TODO: re-enable as soon as the refl-cpp macro issue is fixed https://github.com/veselink1/refl-cpp/issues/63
+    option(WARNINGS_AS_ERRORS "Treat compiler warnings as errors" FALSE) #TODO: re-enable as soon as the refl-cpp macro issue is fixed https://github.com/veselink1/refl-cpp/issues/63
 
-  set(MSVC_WARNINGS
-      /W4 # Baseline reasonable warnings
-      /w14242 # 'identifier': conversion from 'type1' to 'type1', possible loss of data
-      /w14254 # 'operator': conversion from 'type1:field_bits' to 'type2:field_bits', possible loss of data
-      /w14263 # 'function': member function does not override any base class virtual member function
-      /w14265 # 'classname': class has virtual functions, but destructor is not virtual instances of this class may not
-              # be destructed correctly
-      /w14287 # 'operator': unsigned/negative constant mismatch
-      /we4289 # nonstandard extension used: 'variable': loop control variable declared in the for-loop is used outside
-              # the for-loop scope
-      /w14296 # 'operator': expression is always 'boolean_value'
-      /w14311 # 'variable': pointer truncation from 'type1' to 'type2'
-      /w14545 # expression before comma evaluates to a function which is missing an argument list
-      /w14546 # function call before comma missing argument list
-      /w14547 # 'operator': operator before comma has no effect; expected operator with side-effect
-      /w14549 # 'operator': operator before comma has no effect; did you intend 'operator'?
-      /w14555 # expression has no effect; expected expression with side- effect
-      /w14619 # pragma warning: there is no warning number 'number'
-      /w14640 # Enable warning on thread un-safe static member initialization
-      /w14826 # Conversion from 'type1' to 'type_2' is sign-extended. This may cause unexpected runtime behavior.
-      /w14905 # wide string literal cast to 'LPSTR'
-      /w14906 # string literal cast to 'LPWSTR'
-      /w14928 # illegal copy-initialization; more than one user-defined conversion has been implicitly applied
-      /permissive- # standards conformance mode for MSVC compiler.
-  )
+    set(MSVC_WARNINGS
+            /W4 # Baseline reasonable warnings
+            /w14242 # 'identifier': conversion from 'type1' to 'type1', possible loss of data
+            /w14254 # 'operator': conversion from 'type1:field_bits' to 'type2:field_bits', possible loss of data
+            /w14263 # 'function': member function does not override any base class virtual member function
+            /w14265 # 'classname': class has virtual functions, but destructor is not virtual instances of this class may not
+            # be destructed correctly
+            /w14287 # 'operator': unsigned/negative constant mismatch
+            /we4289 # nonstandard extension used: 'variable': loop control variable declared in the for-loop is used outside
+            # the for-loop scope
+            /w14296 # 'operator': expression is always 'boolean_value'
+            /w14311 # 'variable': pointer truncation from 'type1' to 'type2'
+            /w14545 # expression before comma evaluates to a function which is missing an argument list
+            /w14546 # function call before comma missing argument list
+            /w14547 # 'operator': operator before comma has no effect; expected operator with side-effect
+            /w14549 # 'operator': operator before comma has no effect; did you intend 'operator'?
+            /w14555 # expression has no effect; expected expression with side- effect
+            /w14619 # pragma warning: there is no warning number 'number'
+            /w14640 # Enable warning on thread un-safe static member initialization
+            /w14826 # Conversion from 'type1' to 'type_2' is sign-extended. This may cause unexpected runtime behavior.
+            /w14905 # wide string literal cast to 'LPSTR'
+            /w14906 # string literal cast to 'LPWSTR'
+            /w14928 # illegal copy-initialization; more than one user-defined conversion has been implicitly applied
+            /permissive- # standards conformance mode for MSVC compiler.
+    )
 
-  set(CLANG_WARNINGS
-      # -Werror # avoid warnings since they are often indicative of immature API and/or potential sources of bugs # TODO: enable
-      -Wall
-      -Wextra # reasonable and standard
-      -Wshadow # warn the user if a variable declaration shadows one from a parent context
-      -Wnon-virtual-dtor # warn the user if a class with virtual functions has a non-virtual destructor. This helps
-                         # catch hard to track down memory errors
-      -Wold-style-cast # warn for c-style casts
-      -Wcast-align # warn for potential performance problem casts
-      -Wunused # warn on anything being unused
-      -Woverloaded-virtual # warn if you overload (not override) a virtual function
-      -Wpedantic # warn if non-standard C++ is used
-      -Wconversion # warn on type conversions that may lose data
-      -Wsign-conversion # warn on sign conversions
-      -Wnull-dereference # warn if a null dereference is detected
-      -Wdouble-promotion # warn if float is implicit promoted to double
-      -Wformat=2 # warn on security issues around functions that format output (ie printf)
-      -Wno-unknown-pragmas # ignore IDE, GCC/CLANG specific pragmas
-      -Wimplicit-fallthrough # Warns when case statements fall-through.
-  )
+    set(CLANG_WARNINGS
+            # -Werror # avoid warnings since they are often indicative of immature API and/or potential sources of bugs # TODO: enable
+            -Wall
+            -Wextra # reasonable and standard
+            -Wshadow # warn the user if a variable declaration shadows one from a parent context
+            -Wnon-virtual-dtor # warn the user if a class with virtual functions has a non-virtual destructor. This helps
+            # catch hard to track down memory errors
+            -Wold-style-cast # warn for c-style casts
+            -Wcast-align # warn for potential performance problem casts
+            -Wunused # warn on anything being unused
+            -Woverloaded-virtual # warn if you overload (not override) a virtual function
+            -Wpedantic # warn if non-standard C++ is used
+            -Wconversion # warn on type conversions that may lose data
+            -Wsign-conversion # warn on sign conversions
+            -Wnull-dereference # warn if a null dereference is detected
+            -Wdouble-promotion # warn if float is implicit promoted to double
+            -Wformat=2 # warn on security issues around functions that format output (ie printf)
+            -Wno-unknown-pragmas # ignore IDE, GCC/CLANG specific pragmas
+            -Wimplicit-fallthrough # Warns when case statements fall-through.
+    )
 
-  if(WARNINGS_AS_ERRORS)
-    set(CLANG_WARNINGS ${CLANG_WARNINGS} -Werror)
-    set(MSVC_WARNINGS ${MSVC_WARNINGS} /WX)
-  endif()
+    if (WARNINGS_AS_ERRORS)
+        set(CLANG_WARNINGS ${CLANG_WARNINGS} -Werror)
+        set(MSVC_WARNINGS ${MSVC_WARNINGS} /WX)
+    endif ()
 
-  set(GCC_WARNINGS
-          ${CLANG_WARNINGS}
-          -Wno-dangling-reference # TODO: remove this once the fmt dangling reference bug is fixed
-          -Wmisleading-indentation # warn if indentation implies blocks where blocks do not exist
-          -Wduplicated-cond # warn if if / else chain has duplicated conditions
-          -Wduplicated-branches # warn if if / else branches have duplicated code
-          -Wlogical-op # warn about logical operations being used where bitwise were probably wanted
-          -Wuseless-cast # warn if you perform a cast to the same type
-          -Wno-interference-size # suppress ABI compatibility warnings for hardware inferred size
-          -Wno-maybe-uninitialized # false positives if asan is enabled: https://gcc.gnu.org/bugzilla//show_bug.cgi?id=1056h6
-          -fconcepts-diagnostics-depth=3
-          )
+    set(GCC_WARNINGS
+            ${CLANG_WARNINGS}
+            -Wno-dangling-reference # TODO: remove this once the fmt dangling reference bug is fixed
+            -Wmisleading-indentation # warn if indentation implies blocks where blocks do not exist
+            -Wduplicated-cond # warn if if / else chain has duplicated conditions
+            -Wduplicated-branches # warn if if / else branches have duplicated code
+            -Wlogical-op # warn about logical operations being used where bitwise were probably wanted
+            -Wuseless-cast # warn if you perform a cast to the same type
+            -Wno-interference-size # suppress ABI compatibility warnings for hardware inferred size
+            -Wno-maybe-uninitialized # false positives if asan is enabled: https://gcc.gnu.org/bugzilla//show_bug.cgi?id=1056h6
+            -fconcepts-diagnostics-depth=3
+    )
 
-  if(MSVC)
-    set(PROJECT_WARNINGS ${MSVC_WARNINGS})
-  elseif(CMAKE_CXX_COMPILER_ID MATCHES ".*Clang")
-    set(PROJECT_WARNINGS ${CLANG_WARNINGS})
-  elseif(CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
-    set(PROJECT_WARNINGS ${GCC_WARNINGS})
-  else()
-    message(AUTHOR_WARNING "No compiler warnings set for '${CMAKE_CXX_COMPILER_ID}' compiler.")
-  endif()
+    if (MSVC)
+        set(PROJECT_WARNINGS ${MSVC_WARNINGS})
+    elseif (CMAKE_CXX_COMPILER_ID MATCHES ".*Clang")
+        set(PROJECT_WARNINGS ${CLANG_WARNINGS})
+    elseif (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
+        set(PROJECT_WARNINGS ${GCC_WARNINGS})
+    else ()
+        message(AUTHOR_WARNING "No compiler warnings set for '${CMAKE_CXX_COMPILER_ID}' compiler.")
+    endif ()
 
-  target_compile_options(${project_name} INTERFACE ${PROJECT_WARNINGS})
+    target_compile_options(${project_name} INTERFACE ${PROJECT_WARNINGS})
 
 endfunction()