inplace span support (#84)

Use `inplace' to explicitly specify calculations on the input data.

- Add `inplace_span_t<T>` class;
- Add `inplace` function;
- Add some inplace functions;
- Add inplace version of `solve` for `FftPlanC`;
- Code refactoring;
- Remove deprecated;
- Enable cache for fftw (thread-safe problems);

Author: vitalsong

CMakeLists.txt

@@ -71,7 +71,6 @@ if (DSPLIB_FFT_BACKEND STREQUAL "fftw")
     endif()
     set(FFT_INCLUDES ${FFTW_INCLUDES})
     set(DSPLIB_EXCLUDE_FFT ON)
-    set(DSPLIB_FFT_CACHE_SIZE 0) # FFTW has its own cache
 elseif (DSPLIB_FFT_BACKEND STREQUAL "ne10")
     find_package(NE10 REQUIRED)
     if (DSPLIB_USE_FLOAT32)

include/dsplib/czt.h

@@ -18,10 +18,6 @@ class CztPlan : public FftPlanC
 
     [[nodiscard]] int size() const noexcept final;
 
-    arr_cmplx operator()(span_t<cmplx_t> x) const {
-        return this->solve(x);
-    }
-
 private:
     std::shared_ptr<CztPlanImpl> _d;
 };

include/dsplib/fft.h

@@ -25,6 +25,14 @@ class FftPlanC
         r = this->solve(x);
     }
 
+    //inplace FFT implementation
+    virtual void solve(inplace_cmplx inp) const {
+        //default non optimal implementation with temp array
+        auto x = inp.get();
+        const auto y = this->solve(x);
+        x.assign(y);
+    }
+
     [[nodiscard]] virtual int size() const noexcept = 0;
 };
 

include/dsplib/fir.h

@@ -18,20 +18,22 @@ class FirFilter
       , _d(h.size() - 1) {
     }
 
-    //main processing
-    base_array<T> process(const base_array<T>& s) {
-        if (_h.size() == 1) {
-            return s * _h[0];
-        }
-
-        auto x = _d | s;
-        auto r = FirFilter::conv(x, _h);
-        int nd = _h.size() - 1;
-        int nx = x.size();
-        _d = x.slice((nx - nd), nx);
+    base_array<T> process(const span_t<T>& s) {
+        base_array<T> r(s);
+        this->process(inplace(r));
         return r;
     }
 
+    void process(inplace_span_t<T> si) {
+        auto s = si.get();
+        auto x = concatenate(_d, s);
+        FirFilter::conv(x, _h);
+        s.assign(x.slice(0, s.size()));
+        const int nd = _d.size();
+        const int nx = x.size();
+        _d.slice(0, nd) = x.slice((nx - nd), nx);
+    }
+
     //current impulse response
     [[nodiscard]] const base_array<T>& coeffs() const {
         return _h;
@@ -41,14 +43,13 @@ class FirFilter
         return _h;
     }
 
-    //convolution operation
-    static base_array<T> conv(const base_array<T>& x, const base_array<T>& h);
-
-    base_array<T> operator()(const base_array<T>& x) {
+    base_array<T> operator()(const span_t<T>& x) {
         return this->process(x);
     }
 
 private:
+    static void conv(mut_span_t<T> x, span_t<T> h);
+
     base_array<T> _h;   ///< impulse response
     base_array<T> _d;   ///< filter delay
 };

include/dsplib/iterator.h

@@ -55,14 +55,6 @@ struct SliceIterator
         return lhs.ptr_ != rhs.ptr_;
     };
 
-    friend bool operator>(const SliceIterator& lhs, const SliceIterator& rhs) noexcept {
-        return lhs.ptr_ > rhs.ptr_;
-    };
-
-    friend bool operator<(const SliceIterator& lhs, const SliceIterator& rhs) noexcept {
-        return lhs.ptr_ < rhs.ptr_;
-    };
-
 private:
     pointer ptr_;
     int step_{1};

include/dsplib/math.h

@@ -7,6 +7,8 @@
 
 namespace dsplib {
 
+//TODO: mark noexcept
+
 //exponential
 arr_real exp(const arr_real& arr);
 real_t exp(real_t v);
@@ -53,20 +55,23 @@ int argmin(span_real arr);
 int argmin(span_cmplx arr);
 
 //absolute value and complex magnitude
-arr_real abs(const arr_real& arr);
-real_t abs(real_t v);
-arr_real abs(const arr_cmplx& arr);
-real_t abs(cmplx_t v);
+arr_real abs(const arr_real& arr) noexcept;
+real_t abs(real_t v) noexcept;
+arr_real abs(const arr_cmplx& arr) noexcept;
+real_t abs(cmplx_t v) noexcept;
+void abs(inplace_real arr) noexcept;
 
 //phase angle in the interval [-pi, pi] for each element of a complex array z
 arr_real angle(const arr_cmplx& arr);
 real_t angle(cmplx_t v);
 
 //round
-real_t round(const real_t& x);
-cmplx_t round(const cmplx_t& x);
-arr_real round(const arr_real& arr);
-arr_cmplx round(const arr_cmplx& arr);
+real_t round(const real_t& x) noexcept;
+cmplx_t round(const cmplx_t& x) noexcept;
+arr_real round(const arr_real& arr) noexcept;
+arr_cmplx round(const arr_cmplx& arr) noexcept;
+void round(inplace_real arr) noexcept;
+void round(inplace_cmplx arr) noexcept;
 
 //array sum
 real_t sum(const arr_real& arr);
@@ -113,7 +118,9 @@ arr_real imag(const arr_cmplx& x);
 real_t imag(const cmplx_t& x);
 
 //complex pairing
-arr_cmplx conj(const arr_cmplx& x);
+arr_cmplx conj(const arr_cmplx& x) noexcept;
+
+void conj(inplace_cmplx x) noexcept;
 
 constexpr cmplx_t conj(const cmplx_t& x) noexcept {
     return x.conj();
@@ -181,8 +188,9 @@ arr_cmplx power(const arr_cmplx& x, int n);
 
 //square root (only positive values)
 //TODO: add complex result for negative or complex input
-real_t sqrt(real_t x);
-arr_real sqrt(const arr_real& x);
+real_t sqrt(real_t x) noexcept;
+arr_real sqrt(const arr_real& arr) noexcept;
+void sqrt(inplace_real arr) noexcept;
 
 //array log
 arr_real log(const arr_real& arr);
@@ -210,7 +218,7 @@ arr_real upsample(const arr_real& arr, int n, int phase = 0);
 arr_cmplx upsample(const arr_cmplx& arr, int n, int phase = 0);
 
 //abs(x)^2
-arr_real abs2(const arr_cmplx& x);
+arr_real abs2(const arr_cmplx& x) noexcept;
 
 constexpr real_t abs2(const cmplx_t& x) noexcept {
     return x.abs2();
@@ -269,16 +277,20 @@ inline cmplx_t sign(const cmplx_t& x) noexcept {
 
 //----------------------------------------------------------------------------------------
 //convert power <-> decibels: db = 10 * log10(pow)
-real_t pow2db(real_t v);
-arr_real pow2db(const arr_real& v);
-real_t db2pow(real_t v);
-arr_real db2pow(const arr_real& v);
+real_t pow2db(real_t v) noexcept;
+arr_real pow2db(const arr_real& arr) noexcept;
+real_t db2pow(real_t v) noexcept;
+arr_real db2pow(const arr_real& arr) noexcept;
+void pow2db(inplace_real arr) noexcept;
+void db2pow(inplace_real arr) noexcept;
 
 //convert magnitude <-> decibels: db = 20 * log10(mag)
-real_t mag2db(real_t v);
-arr_real mag2db(const arr_real& v);
-real_t db2mag(real_t v);
-arr_real db2mag(const arr_real& v);
+real_t mag2db(real_t v) noexcept;
+arr_real mag2db(const arr_real& arr) noexcept;
+real_t db2mag(real_t v) noexcept;
+arr_real db2mag(const arr_real& arr) noexcept;
+void mag2db(inplace_real arr) noexcept;
+void db2mag(inplace_real arr) noexcept;
 
 //----------------------------------------------------------------------------------------
 //check that the number is prime

include/dsplib/slice.h

@@ -120,17 +120,17 @@ class mut_slice_t
     }
 
     mut_slice_t& operator=(const slice_t<T>& rhs) {
-        this->copy(*this, rhs);
+        this->assign(rhs);
         return *this;
     }
 
     mut_slice_t& operator=(const mut_slice_t<T>& rhs) {
-        *this = slice_t<T>(rhs);
+        this->assign(slice_t<T>(rhs));
         return *this;
     }
 
     mut_slice_t& operator=(const base_array<T>& rhs) {
-        DSPLIB_ASSERT(!is_same_memory(rhs.slice(0, rhs.size()), *this), "Assigned array to same slice");
+        DSPLIB_ASSERT(!is_same_memory(rhs.slice(0, rhs.size())), "Assigned array to same slice");
         *this = rhs.slice(0, rhs.size());
         return *this;
     }
@@ -165,19 +165,42 @@ class mut_slice_t
         return current;
     }
 
+    //TODO: replace to `copy` or `to_arr` function
     base_array<T> operator*() const noexcept {
         return base_array<T>(*this);
     }
 
-    static bool is_same_memory(slice_t<T> r1, slice_t<T> r2) noexcept {
-        if (r1.empty() || r2.empty()) {
-            return false;
+    void assign(slice_t<T> rhs) {
+        DSPLIB_ASSERT(size() == rhs.size(), "Slices size must be equal");
+        const int count = size();
+
+        //empty slice assign
+        if (count == 0) {
+            return;
         }
-        auto start1 = r1.begin();
-        auto end1 = r1.end();
-        auto start2 = r2.begin();
-        auto end2 = r2.end();
-        return (start1 < end2) && (start2 < end1);
+
+        //simple block copy/move (optimization)
+        const bool is_same = is_same_memory(rhs);
+
+        //check all slices is span
+        if ((stride() == 1) && (rhs.stride() == 1)) {
+            const auto* src = rhs.data_;
+            auto* dst = data_;
+            if (!is_same) {
+                std::memcpy(dst, src, count * sizeof(*src));
+            } else {
+                std::memmove(dst, src, count * sizeof(*src));
+            }
+            return;
+        }
+
+        //same array, specific indexing
+        if (is_same) {
+            *this = base_array<T>(rhs);
+            return;
+        }
+
+        std::copy(rhs.begin(), rhs.end(), begin());
     }
 
     static mut_slice_t make_slice(T* data, int size, int i1, int i2, int step) {
@@ -202,38 +225,6 @@ class mut_slice_t
         return mut_slice_t(data + i1, step, count);
     }
 
-    static void copy(mut_slice_t<T> lhs, slice_t<T> rhs) {
-        DSPLIB_ASSERT(lhs.size() == rhs.size(), "Slices size must be equal");
-        const int count = lhs.size();
-
-        //empty slice assign
-        if (count == 0) {
-            return;
-        }
-
-        //simple block copy/move (optimization)
-        const bool is_same = is_same_memory(rhs, slice_t(lhs));
-        if ((lhs.stride() == 1) && (rhs.stride() == 1)) {
-            const auto* src = &(*rhs.begin());
-            auto* dst = &(*lhs.begin());
-            if (!is_same) {
-                std::memcpy(dst, src, count * sizeof(*src));
-            } else {
-                //overlapped
-                std::memmove(dst, src, count * sizeof(*src));
-            }
-            return;
-        }
-
-        //same array, specific indexing
-        if (is_same) {
-            lhs = base_array<T>(rhs);
-            return;
-        }
-
-        std::copy(rhs.begin(), rhs.end(), lhs.begin());
-    }
-
 protected:
     explicit mut_slice_t(T* data, int stride, int count)
       : data_{data}
@@ -242,6 +233,25 @@ class mut_slice_t
         DSPLIB_ASSERT(count >= 0, "Count of elements must be positive");
     }
 
+    bool is_same_memory(slice_t<T> rhs) noexcept {
+        if (empty() || rhs.empty()) {
+            return false;
+        }
+        auto start1 = rhs.data_;
+        auto end1 = start1 + (rhs.count_ * rhs.stride_);
+        if (start1 > end1) {
+            std::swap(start1, end1);
+        }
+
+        auto start2 = data_;
+        auto end2 = start2 + (count_ * stride_);
+        if (start2 > end2) {
+            std::swap(start2, end2);
+        }
+
+        return (start1 < end2) && (start2 < end1);
+    }
+
     T* data_{nullptr};
     int stride_{0};
     int count_{0};