Feature/real fft

.github/workflows/linux.yml

@@ -22,7 +22,7 @@ jobs:
       - uses: actions/checkout@v2
 
       - name: Configure CMake
-        run: cmake -B ${{github.workspace}}/build-${{matrix.compiler}} -DDSPLIB_BUILD_TESTS=ON -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DDSPLIB_ASAN_ENABLED=True
+        run: cmake -B ${{github.workspace}}/build-${{matrix.compiler}} -DDSPLIB_BUILD_TESTS=ON -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DDSPLIB_ASAN_ENABLED=ON
         env:
           CXX: ${{matrix.compiler}}
 

README.md

@@ -67,6 +67,15 @@ arr_cmplx y1 = fft(x);  //500
 arr_cmplx y2 = fft(x, 1024); //1024
 ```
 
+### FFT performance (double complex)
+`Intel Core i7 (2.6 GHz)`
+
+| N    | dsplib   | fftw3    | kissfft  |
+| ---- | -------- | -------- | -------  |
+| 2048 | 15651 ns | 4357 ns  | 19846 ns |
+| 4096 | 34818 ns | 11726 ns | 36083 ns |
+| 8192 | 81687 ns | 29223 ns | 100624 ns |
+
 ### Inverse Fast Fourier Transform:
 ```cpp
 arr_cmplx x = 1i * zeros(512);

benchs/CMakeLists.txt

@@ -5,6 +5,7 @@ set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 include(${CMAKE_SOURCE_DIR}/cmake/CPM.cmake)
+include(${CMAKE_SOURCE_DIR}/cmake/FindFFTW3.cmake)
 
 cpmaddpackage(NAME benchmark
   GIT_REPOSITORY "https://github.com/google/benchmark.git"
@@ -17,11 +18,27 @@ cpmaddpackage(NAME benchmark
   "BENCHMARK_INSTALL_DOCS OFF"
 )
 
+cpmaddpackage(NAME kissfft
+  GIT_REPOSITORY "https://github.com/mborgerding/kissfft.git"
+  GIT_TAG "131.1.0"
+  OPTIONS
+  "KISSFFT_DATATYPE double"
+  "KISSFFT_TEST OFF"
+  "KISSFFT_STATIC ON"
+)
+
 set(SOURCES
     main.cpp
+    fft.cpp
     adaptive.cpp
 )
 
 add_executable(${PROJECT_NAME} ${SOURCES})
 target_include_directories(${PROJECT_NAME} PUBLIC ${CMAKE_CURRENT_LIST_DIR})
-target_link_libraries(${PROJECT_NAME} dsplib benchmark)
+target_link_libraries(${PROJECT_NAME} dsplib benchmark kissfft)
+
+# CheckFFTW3()
+# if (FFTW_DOUBLE_LIB_FOUND)
+#   target_link_libraries(${PROJECT_NAME} FFT:FFTW3_DOUBLE)
+#   target_compile_definitions(${PROJECT_NAME} PUBLIC FFTW3_SUPPORT)
+# endif()

include/dsplib.h

@@ -2,7 +2,6 @@
 
 #include <dsplib/array.h>
 #include <dsplib/fft.h>
-#include <dsplib/ifft.h>
 #include <dsplib/hilbert.h>
 #include <dsplib/fir.h>
 #include <dsplib/math.h>

include/dsplib/czt.h

@@ -12,7 +12,7 @@ class CztPlan
 public:
     explicit CztPlan(int n, int m, cmplx_t w, cmplx_t a = 1);
     arr_cmplx operator()(const arr_cmplx& x) const;
-    arr_cmplx solve(const arr_cmplx& x) const;
+    [[nodiscard]] arr_cmplx solve(const arr_cmplx& x) const;
 
 private:
     std::shared_ptr<CztPlanImpl> _d;

include/dsplib/fft.h

@@ -7,21 +7,38 @@ namespace dsplib {
 
 class FftPlanImpl;
 
+//Fast Fourier Transform calculation plan
 class FftPlan
 {
 public:
-    FftPlan(int n);
-    arr_cmplx operator()(const arr_cmplx& x) const;
+    explicit FftPlan(int n);
+    [[nodiscard]] arr_cmplx operator()(const arr_cmplx& x) const {
+        return this->solve(x);
+    }
     [[nodiscard]] arr_cmplx solve(const arr_cmplx& x) const;
     [[nodiscard]] int size() const noexcept;
 
 private:
     std::shared_ptr<FftPlanImpl> _d;
 };
 
-using fft_plan [[deprecated]] = FftPlan;
+//Inverse Fast Fourier Transform calculation plan
+class IFftPlan : public FftPlan
+{
+public:
+    explicit IFftPlan(int n)
+      : FftPlan(n) {
+    }
 
-//TODO: rfft optimization
+    [[nodiscard]] arr_cmplx operator()(const arr_cmplx& x) const {
+        return this->solve(x);
+    }
+    [[nodiscard]] arr_cmplx solve(const arr_cmplx& x) const;
+    [[nodiscard]] int size() const noexcept;
+};
+
+using fft_plan [[deprecated]] = FftPlan;
+using ifft_plan [[deprecated]] = IFftPlan;
 
 /*!
  * \brief Fast Fourier Transform (complex)
@@ -36,4 +53,15 @@ arr_cmplx fft(const arr_cmplx& x);
 // if x.size() is greater than n, then x is truncated to length n
 arr_cmplx fft(const arr_cmplx& x, int n);
 
+//Fast Fourier Transform (real)
+arr_cmplx fft(const arr_real& x);
+
+/*!
+ * \brief Inverse fourier transform
+ * \details IFFT for complex signal
+ * \param x Input array [N]
+ * \return Result array [N]
+ */
+arr_cmplx ifft(const arr_cmplx& x);
+
 }   // namespace dsplib

include/dsplib/math.h

@@ -89,7 +89,7 @@ arr_real imag(const arr_cmplx& x);
 real_t imag(cmplx_t x);
 
 //complex pairing
-arr_cmplx conj(const arr_cmplx& x);
+arr_cmplx conj(arr_cmplx x);
 
 constexpr cmplx_t conj(cmplx_t x) {
     return x.conj();
@@ -110,8 +110,8 @@ arr_cmplx complex(const arr_real& re, const arr_real& im);
 int nextpow2(int m);
 
 //array pow
-arr_real pow2(const arr_real& arr);
-arr_cmplx pow2(const arr_cmplx& arr);
+arr_real pow2(arr_real arr);
+arr_cmplx pow2(arr_cmplx arr);
 
 constexpr real_t pow2(real_t x) {
     return x * x;

include/dsplib/utils.h

@@ -81,8 +81,8 @@ arr_real repelem(const arr_real& x, int n);
 arr_cmplx repelem(const arr_cmplx& x, int n);
 
 //flip order of elements
-arr_real flip(const arr_real& x);
-arr_cmplx flip(const arr_cmplx& x);
+arr_real flip(arr_real x);
+arr_cmplx flip(arr_cmplx x);
 
 enum class dtype
 {
@@ -201,4 +201,8 @@ real_t peakloc(const arr_cmplx& x, int idx, bool cyclic = true);
 int finddelay(const dsplib::arr_real& x1, const dsplib::arr_real& x2);
 int finddelay(const dsplib::arr_cmplx& x1, const dsplib::arr_cmplx& x2);
 
+//Shift array circularly
+arr_real circshift(const arr_real& x, int shift);
+arr_cmplx circshift(const arr_cmplx& x, int shift);
+
 }   // namespace dsplib

lib/czt.cpp

@@ -1,6 +1,5 @@
 #include <dsplib/czt.h>
 #include <dsplib/fft.h>
-#include <dsplib/ifft.h>
 #include <dsplib/math.h>
 #include <dsplib/utils.h>
 
@@ -23,11 +22,16 @@ class CztPlanImpl
         _rp = chirp.slice(_n - 1, _m + _n - 1);
     }
 
-    arr_cmplx solve(const arr_cmplx& x) const {
-        auto xp = (x * _cp) | zeros(_ich.size() - x.size());
-        auto r = ifft(fft(xp) * _ich);
-        arr_cmplx tr = r.slice(_n - 1, _m + _n - 1);
-        return tr * _rp;
+    [[nodiscard]] arr_cmplx solve(arr_cmplx x) const noexcept {
+        assert(x.size() == _m);
+        x *= _cp;
+        auto z = fft(x, _ich.size());
+        z *= _ich;
+        const auto r = ifft(z);
+        for (int i = 0; i < _m; ++i) {
+            x[i] = r[_n - 1 + i] * _rp[i];
+        }
+        return x;
     }
 
     int _n;