Hide IFFT to factory with caching

Author: vitalsong

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.10)
-project(dsplib LANGUAGES CXX VERSION 0.54.8)
+project(dsplib LANGUAGES CXX VERSION 0.54.9)
 
 set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)

include/dsplib/ifft.h

@@ -1,36 +1,48 @@
 #pragma once
 
-#include <dsplib/types.h>
 #include <dsplib/array.h>
-#include <dsplib/fft.h>
 #include <memory>
 
 namespace dsplib {
 
-class IfftPlan
+class BaseIfftPlanC
 {
 public:
-    IfftPlan(int n);
+    virtual ~BaseIfftPlanC() = default;
+    [[nodiscard]] virtual arr_cmplx solve(const arr_cmplx& x) const = 0;
+    [[nodiscard]] virtual int size() const noexcept = 0;
+};
+
+class BaseIfftPlanR
+{
+public:
+    virtual ~BaseIfftPlanR() = default;
+    [[nodiscard]] virtual arr_real solve(const arr_cmplx& x) const = 0;
+    [[nodiscard]] virtual int size() const noexcept = 0;
+};
+
+class IfftPlan : public BaseIfftPlanC
+{
+public:
+    explicit IfftPlan(int n);
     arr_cmplx operator()(const arr_cmplx& x) const;
-    [[nodiscard]] arr_cmplx solve(const arr_cmplx& x) const;
-    [[nodiscard]] int size() const noexcept;
+    [[nodiscard]] arr_cmplx solve(const arr_cmplx& x) const final;
+    [[nodiscard]] int size() const noexcept final;
 
 private:
-    std::shared_ptr<FftPlan> _d;
+    std::shared_ptr<BaseIfftPlanC> ifft_;
 };
 
-class IfftPlanR
+class IfftPlanR : public BaseIfftPlanR
 {
 public:
-    IfftPlanR(int n);
+    explicit IfftPlanR(int n);
     arr_real operator()(const arr_cmplx& x) const;
-    [[nodiscard]] arr_real solve(const arr_cmplx& x) const;
-    [[nodiscard]] int size() const noexcept;
+    [[nodiscard]] arr_real solve(const arr_cmplx& x) const final;
+    [[nodiscard]] int size() const noexcept final;
 
 private:
-    const int _n;
-    std::shared_ptr<FftPlan> _d;
-    const std::vector<cmplx_t> _w;   //exp(1i * 2 * pi / n)
+    std::shared_ptr<BaseIfftPlanR> ifft_;
 };
 
 /*!

lib/fft/cmplx-ifft.h

@@ -0,0 +1,39 @@
+#pragma once
+
+#include <dsplib/ifft.h>
+
+#include "fft/factory.h"
+
+namespace dsplib {
+
+class CmplxIfftPlan : public BaseIfftPlanC
+{
+public:
+    explicit CmplxIfftPlan(int n)
+      : fft_{create_fft_plan(n)} {
+    }
+
+    arr_cmplx solve(const arr_cmplx& x) const final {
+        const real_t m = real_t(1) / x.size();
+        arr_cmplx y = x * m;
+        _inplace_conj(y);
+        y = fft_->solve(y);
+        _inplace_conj(y);
+        return y;
+    }
+
+    int size() const noexcept final {
+        return fft_->size();
+    }
+
+private:
+    static void _inplace_conj(arr_cmplx& x) {
+        for (auto& v : x) {
+            v.im = -v.im;
+        }
+    }
+
+    std::shared_ptr<BaseFftPlanC> fft_;
+};
+
+}   // namespace dsplib

lib/fft/factory.cpp

@@ -1,12 +1,15 @@
 #include <dsplib/fft.h>
+#include <dsplib/ifft.h>
 #include <dsplib/math.h>
 #include <dsplib/utils.h>
 
+#include "fft/cmplx-ifft.h"
 #include "fft/fact-fft.h"
 #include "fft/primes-fft.h"
 #include "fft/factory.h"
 #include "fft/pow2-fft.h"
 #include "fft/real-fft.h"
+#include "fft/real-ifft.h"
 #include "fft/small-fft.h"
 #include "lru-cache.h"
 
@@ -41,6 +44,14 @@ std::shared_ptr<BaseFftPlanR> _get_rfft_plan(int n) {
     return std::make_shared<FactorFFTPlanR>(n);
 }
 
+std::shared_ptr<BaseIfftPlanC> _get_ifft_plan(int n) {
+    return std::make_shared<CmplxIfftPlan>(n);
+}
+
+std::shared_ptr<BaseIfftPlanR> _get_irfft_plan(int n) {
+    return std::make_shared<RealIfftPlan>(n);
+}
+
 }   // namespace
 
 //-------------------------------------------------------------------------------------------------
@@ -74,4 +85,24 @@ std::shared_ptr<BaseFftPlanR> create_rfft_plan(int n) {
     return cache.get(n);
 }
 
+std::shared_ptr<BaseIfftPlanC> create_ifft_plan(int n) {
+    thread_local LRUCache<int, std::shared_ptr<BaseIfftPlanC>> cache{FFT_CACHE_SIZE};
+    if (!cache.exists(n)) {
+        auto plan = _get_ifft_plan(n);
+        cache.put(n, plan);
+        return plan;
+    }
+    return cache.get(n);
+}
+
+std::shared_ptr<BaseIfftPlanR> create_irfft_plan(int n) {
+    thread_local LRUCache<int, std::shared_ptr<BaseIfftPlanR>> cache{FFT_CACHE_SIZE};
+    if (!cache.exists(n)) {
+        auto plan = _get_irfft_plan(n);
+        cache.put(n, plan);
+        return plan;
+    }
+    return cache.get(n);
+}
+
 }   // namespace dsplib

lib/fft/factory.h

@@ -1,13 +1,18 @@
 #pragma once
 
+#include <dsplib/ifft.h>
 #include <dsplib/fft.h>
 
 namespace dsplib {
 
 //instance or get cached fft plan
 
+//FFT
 std::shared_ptr<BaseFftPlanC> create_fft_plan(int n);
-
 std::shared_ptr<BaseFftPlanR> create_rfft_plan(int n);
 
+//IFFT
+std::shared_ptr<BaseIfftPlanC> create_ifft_plan(int n);
+std::shared_ptr<BaseIfftPlanR> create_irfft_plan(int n);
+
 }   // namespace dsplib

lib/fft/fft.cpp

@@ -4,8 +4,6 @@
 
 #include "fft/factory.h"
 
-#include <cassert>
-
 namespace dsplib {
 
 //-------------------------------------------------------------------------------------------------

lib/fft/ifft.cpp

@@ -2,95 +2,42 @@
 #include <dsplib/fft.h>
 #include <dsplib/math.h>
 
-namespace dsplib {
-
-namespace {
+#include "fft/factory.h"
 
-void _inplace_conj(arr_cmplx& x) {
-    for (auto& v : x) {
-        v.im = -v.im;
-    }
-}
-
-std::vector<cmplx_t> _irfft_coeffs(int n) noexcept {
-    assert(n % 4 == 0);
-    const int n4 = n / 4;
-    const int n2 = n / 2;
-    std::vector<cmplx_t> res(n / 2);
-    res[0] = {1, 0};
-    res[n4] = {0, 1};
-    //use only first n/4 samples
-    for (int i = 1; i < n4; ++i) {
-        const auto v = std::cos(2 * pi * i / n);
-        res[i].re = v;
-        res[n2 - i].re = -v;
-        res[n4 + i].im = v;
-        res[n4 - i].im = v;
-    }
-    return res;
-}
-
-}   // namespace
+namespace dsplib {
 
 //---------------------------------------------------------------------------------
 IfftPlan::IfftPlan(int n)
-  : _d{std::make_shared<FftPlan>(n)} {
+  : ifft_{create_ifft_plan(n)} {
 }
 
 arr_cmplx IfftPlan::operator()(const arr_cmplx& x) const {
     return this->solve(x);
 }
 
 arr_cmplx IfftPlan::solve(const arr_cmplx& x) const {
-    const real_t m = real_t(1) / x.size();
-    arr_cmplx y = x * m;
-    _inplace_conj(y);
-    y = _d->solve(y);
-    _inplace_conj(y);
-    return y;
+    return ifft_->solve(x);
 }
 
 int IfftPlan::size() const noexcept {
-    return _d->size();
+    return ifft_->size();
 }
 
 //---------------------------------------------------------------------------------
 IfftPlanR::IfftPlanR(int n)
-  : _n{n}
-  , _d{std::make_shared<FftPlan>(n / 2)}
-  , _w(_irfft_coeffs(n)) {
-    DSPLIB_ASSERT(n % 2 == 0, "ifft size must be even");
+  : ifft_{create_irfft_plan(n)} {
 }
 
 arr_real IfftPlanR::operator()(const arr_cmplx& x) const {
     return this->solve(x);
 }
 
 arr_real IfftPlanR::solve(const arr_cmplx& x) const {
-    assert(_n % 2 == 0);
-    DSPLIB_ASSERT((x.size() == _n) || (x.size() == _n / 2 + 1), "input size must be n/2+1 or n");
-
-    const real_t dn = real_t(1) / _n;
-    std::vector<cmplx_t> Z(_n / 2);
-    for (int i = 0; i < _n / 2; ++i) {
-        const auto v = conj(x[_n / 2 - i]);
-        const cmplx_t Xe = (x[i] + v) * dn;
-        const cmplx_t Xo = (x[i] - v) * dn * _w[i];
-        Z[i].re = Xe.re - Xo.im;
-        Z[i].im = -Xe.im - Xo.re;
-    }
-
-    const auto z = _d->solve(Z);
-    std::vector<real_t> r(_n);
-    for (int i = 0; i < _n / 2; ++i) {
-        r[2 * i] = z[i].re;
-        r[2 * i + 1] = -z[i].im;
-    }
-    return r;
+    return ifft_->solve(x);
 }
 
 int IfftPlanR::size() const noexcept {
-    return _n;
+    return ifft_->size();
 }
 
 //---------------------------------------------------------------------------------
@@ -104,7 +51,6 @@ arr_real irfft(const arr_cmplx& x) {
 }
 
 arr_real irfft(const arr_cmplx& x, int n) {
-    //FIXIT: _irfft_coeffs overhead, use cache
     IfftPlanR plan(n);
     return plan(x);
 }

lib/fft/real-fft.cpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <dsplib/fft.h>
+#include <dsplib/utils.h>
+#include <dsplib/math.h>
+
+#include "fft/real-fft.h"
+#include "fft/factory.h"
+
+namespace dsplib {
+
+RealFftPlan::RealFftPlan(int n)
+  : n_{n}
+  , fft_{create_fft_plan(n / 2)}
+  , w_(expj(-2 * pi * arange(n / 2) / n).to_vec()) {   //TODO: optimization (like in real ifft module)
+    DSPLIB_ASSERT(n % 2 == 0, "FFT size must be even");
+}
+
+arr_cmplx RealFftPlan::solve(const arr_real& x) const {
+    using namespace std::complex_literals;
+    DSPLIB_ASSERT(x.size() == n_, "Input size must be equal FFT size");
+    const int n2 = n_ / 2;
+
+    arr_cmplx z(reinterpret_cast<const cmplx_t*>(x.data()), n2);
+    const auto Z = fft_->solve(z * 0.5);
+
+    arr_cmplx res(n_);
+
+    {
+        const auto Xe = Z[0] + conj(Z[0]);
+        const auto Xo = (conj(Z[0]) - Z[0]) * w_[0];
+        res[0].re = Xe.re - Xo.im;
+        res[0].im = Xe.im + Xo.re;
+    }
+
+    for (int i = 1; i < n2; ++i) {
+        const auto Zc = conj(Z[n2 - i]);
+        const auto Xe = Z[i] + Zc;
+        const auto Xo = (Zc - Z[i]) * w_[i];
+        res[i].re = Xe.re - Xo.im;
+        res[i].im = Xe.im + Xo.re;
+        res[n_ - i].re = res[i].re;
+        res[n_ - i].im = -res[i].im;
+    }
+
+    {
+        const auto Xe = Z[0] + conj(Z[0]);
+        const auto Xo = conj(Z[0]) - Z[0];
+        res[n2] = Xe.re + Xo.im;
+    }
+
+    return res;
+}
+
+int RealFftPlan::size() const noexcept {
+    return n_;
+}
+
+}   // namespace dsplib