make FftFilter as template class (#95)

Author: vitalsong

CMakeLists.txt

@@ -24,6 +24,7 @@ set(DSPLIB_SOURCES
     lib/detector.cpp
     lib/findpeaks.cpp
     lib/fir.cpp
+    lib/fft-filter.cpp
     lib/gccphat.cpp
     lib/hilbert.cpp
     lib/math.cpp

include/dsplib/agc.h

@@ -8,6 +8,7 @@ namespace dsplib {
 struct AgcImpl;
 
 //Adaptively adjust gain for constant signal level output
+//TODO: make as template
 class Agc
 {
 public:

include/dsplib/detector.h

@@ -17,6 +17,7 @@ class PreambleDetectorImpl;
 // corr = filter(conj(h) / (nh * rms(h)), 1, x);
 // pagg = filter(ones(1, nh) / nh, 1, abs(x).^2);
 // res = sqrt((abs(corr) .^ 2) ./ pagg);
+//TODO: add real type implementation
 class PreambleDetector
 {
 public:

include/dsplib/fir.h

@@ -2,6 +2,9 @@
 
 #include <dsplib/array.h>
 #include <dsplib/keywords.h>
+#include <dsplib/math.h>
+#include <dsplib/fft.h>
+#include <dsplib/ifft.h>
 
 namespace dsplib {
 
@@ -60,45 +63,37 @@ using FirFilterC = FirFilter<cmplx_t>;
 template<typename U>
 FirFilter(const base_array<U>&) -> FirFilter<U>;
 
+template<typename T>
+class FftFilterImpl;
+
 /*!
  * \brief FFT-based FIR filtering using overlap-add method
- * \details Fast fir implementation for IR len > 200
- * \todo make template (real version is faster)
+ * \details Fast fir implementation for large IR length (usually > 200)
  */
+template<typename T>
 class FftFilter
 {
 public:
-    FftFilter() = default;
-    explicit FftFilter(span_real h);
+    explicit FftFilter(span_t<T> h);
 
-    explicit FftFilter(span_cmplx h);
-
-    //usually in.size() != out.size()
-    arr_real process(span_real x);
-    arr_cmplx process(span_cmplx x);
-
-    //optimal input size for y[nx] = process(x[nx])
-    [[nodiscard]] int block_size() const {
-        return _n;
-    }
+    base_array<T> process(span_t<T> x);
 
-    arr_real operator()(span_real x) {
+    base_array<T> operator()(span_t<T> x) {
         return this->process(x);
     }
 
-    arr_cmplx operator()(span_cmplx x) {
-        return this->process(x);
-    }
+    [[nodiscard]] int block_size() const;
 
 private:
-    arr_cmplx _x;
-    arr_cmplx _h;
-    arr_cmplx _olap;
-    int _nx{0};
-    int _m{0};
-    int _n{0};
+    std::shared_ptr<FftFilterImpl<T>> d_;
 };
 
+using FftFilterR = FftFilter<real_t>;
+using FftFilterC = FftFilter<cmplx_t>;
+
+template<typename U>
+FftFilter(const base_array<U>&) -> FftFilter<U>;
+
 //Type of linear phase FIR filter
 FirType firtype(span_real h);
 

lib/detector.cpp

@@ -89,7 +89,7 @@ class PreambleDetectorImpl
 
     void reset() {
         _pow_flt.process(zeros(frame_len()));
-        _corr_flt.process(zeros(frame_len()));
+        _corr_flt.process(complex(zeros(frame_len())));
         _delay.reset();
     }
 
@@ -98,7 +98,7 @@ class PreambleDetectorImpl
         return flip(h) / (rms(h) * h.size());
     }
 
-    FftFilter _corr_flt;
+    FftFilterC _corr_flt;
     MAFilterR _pow_flt;
     real_t _threshold{1.0};
     CDelay<cmplx_t> _delay;

lib/fft-filter.cpp

@@ -0,0 +1,91 @@
+#include "dsplib/fir.h"
+
+namespace dsplib {
+
+template<typename T>
+class FftFilterImpl
+{
+public:
+    explicit FftFilterImpl(span_t<T> h)
+      : _nfft{int(1) << nextpow2(2 * h.size())}
+      , _m{h.size()}
+      , _n{_nfft - h.size() + 1}
+      , _x(_nfft)
+      , _olap(_m - 1) {
+        assert(_n > _m);
+        _h = fft(conj(h), _nfft);
+    }
+
+    base_array<T> process(span_t<T> x) {
+        const int nr = (x.size() + _nx) / _n * _n;
+        base_array<T> r(nr);
+        auto* pr = r.data();
+        for (const auto& val : x) {
+            _x[_nx] = val;
+            _nx += 1;
+            if (_nx == _n) {
+                base_array<T> ry;
+                if constexpr (std::is_same_v<T, real_t>) {
+                    ry = std::move(irfft(fft(_x) * _h));   //TODO: use n/2+1 multiply
+                } else {
+                    ry = std::move(ifft(fft(_x) * _h));
+                }
+                for (int i = 0; i < _n; i++) {
+                    pr[i] = ry[i];
+                }
+                for (int i = 0; i < (_m - 1); i++) {
+                    pr[i] += _olap[i];
+                    _olap[i] = ry[i + _n];
+                }
+                pr += _n;
+                _nx = 0;
+            }
+        }
+        return r;
+    }
+
+    [[nodiscard]] int block_size() const {
+        return _n;
+    }
+
+private:
+    int _nx{0};
+    int _nfft;
+    int _m{0};
+    int _n{0};
+    base_array<T> _x;
+    arr_cmplx _h;
+    base_array<T> _olap;
+};
+
+template<>
+FftFilter<real_t>::FftFilter(span_t<real_t> h) {
+    d_ = std::make_shared<FftFilterImpl<real_t>>(h);
+}
+
+template<>
+FftFilter<cmplx_t>::FftFilter(span_t<cmplx_t> h) {
+    d_ = std::make_shared<FftFilterImpl<cmplx_t>>(h);
+}
+
+template<>
+base_array<real_t> FftFilter<real_t>::process(span_t<real_t> x) {
+    return d_->process(x);
+}
+
+template<>
+base_array<cmplx_t> FftFilter<cmplx_t>::process(span_t<cmplx_t> x) {
+    return d_->process(x);
+}
+
+template<>
+int FftFilter<real_t>::block_size() const {
+    return d_->block_size();
+}
+
+template<>
+int FftFilter<cmplx_t>::block_size() const {
+    return d_->block_size();
+}
+
+}   // namespace dsplib

lib/fir.cpp

@@ -119,52 +119,6 @@ void FirFilter<cmplx_t>::conv(mut_span_t<cmplx_t> x, span_t<cmplx_t> h) {
     _conv(x.data(), h.data(), h.size(), x.size());
 }
 
-//-------------------------------------------------------------------------------------------------
-FftFilter::FftFilter(span_cmplx h)
-  : _m{h.size()} {
-    const int fft_len = 1L << nextpow2(2 * h.size());
-    _n = fft_len - h.size() + 1;
-    assert(_n > _m);
-    _olap = complex(zeros(_m - 1));
-    _h = fft(conj(h), fft_len);
-    _x = complex(zeros(fft_len));
-}
-
-FftFilter::FftFilter(span_real h)
-  : FftFilter(complex(h)) {
-}
-
-arr_cmplx FftFilter::process(span_cmplx x) {
-    const int nr = (x.size() + _nx) / _n * _n;
-    arr_cmplx r(nr);
-    cmplx_t* pr = r.data();
-    for (const auto& val : x) {
-        _x[_nx] = val;
-        _nx += 1;
-        if (_nx == _n) {
-            const auto ry = ifft(fft(_x) * _h);
-
-            for (int i = 0; i < _n; i++) {
-                pr[i] = ry[i];
-            }
-
-            for (int i = 0; i < (_m - 1); i++) {
-                pr[i] += _olap[i];
-                _olap[i] = ry[i + _n];
-            }
-
-            pr += _n;
-            _nx = 0;
-        }
-    }
-    return r;
-}
-
-arr_real FftFilter::process(span_real x) {
-    //TODO: real optimization
-    return real(process(complex(x)));
-}
-
 //----------------------------------------------------------------------------------------------
 arr_real fir1(int n, real_t wn, FilterType ftype, const arr_real& win) {
     assert(n > 0);