Span API support (#51)

* Span API support

- Remove auto cast `arr_real` to `arr_cmplx`;
- Add `complex(re)` function;
- Add `span` implementation;
- Return `span` from `slice` function for step=1;
- Remove `base_array` reference from `slice`;
- Check same array by pointer range;
- Simplify cast for `cmplx_t` type;
- Add span iterators;
- Add span api for some functions;
- Add span API for `concatenate`. Empty span support;
- Remove base_array::operator(int i);
- Add non-copy fft solve functions;
- Add CMakePresets.json;
- Remove deprecated;
- Rename `slice_t` to `mut_slice_t`;
- Rename `const_slice_t` to `slice_t`;
- Limit the set of types for base_array<T>;
- Rename `span` to `make_span`;
- Add some `DSPLIB_ASSERT`;
- Fix benchs error;
- Code refactoring;
- Fix CI config;
- Update tests;

Author: vitalsong

.github/workflows/android-ndk.yml

@@ -5,7 +5,7 @@ on:
     paths-ignore: ["README.md"]
 
   pull_request:
-    branches: [master]
+    branches: [master, develop]
 
 env:
   BUILD_TYPE: Release

.github/workflows/linux.yml

@@ -5,7 +5,7 @@ on:
     paths-ignore: ["README.md"]
 
   pull_request:
-    branches: [master]
+    branches: [master, develop]
 
 env:
   BUILD_TYPE: Debug

.github/workflows/windows.yml

@@ -5,7 +5,7 @@ on:
     paths-ignore: ["README.md"]
 
   pull_request:
-    branches: [master]
+    branches: [master, develop]
 
 env:
   BUILD_TYPE: Release

CMakeLists.txt

@@ -46,6 +46,7 @@ set(DSPLIB_SOURCES
     lib/ifft.cpp
     lib/czt.cpp
     lib/subband.cpp
+    lib/internal/besseli.cpp
 )
 
 if (DSPLIB_EXCLUDE_FFT)

CMakePresets.json

@@ -0,0 +1,68 @@
+{
+    "version": 8,
+    "cmakeMinimumRequired": {
+        "major": 3,
+        "minor": 18,
+        "patch": 0
+    },
+    "configurePresets": [
+        {
+            "name": "default",
+            "hidden": true,
+            "generator": "Ninja",
+            "binaryDir": "${sourceDir}/build",
+            "cacheVariables": {
+                "CMAKE_BUILD_TYPE": "Debug",
+                "DSPLIB_BUILD_TESTS": "OFF",
+                "DSPLIB_BUILD_BENCHS": "OFF",
+                "DSPLIB_BUILD_EXAMPLES": "OFF",
+                "DSPLIB_USE_FLOAT32": "OFF",
+                "DSPLIB_EXCLUDE_FFT": "OFF",
+                "DSPLIB_ASAN_ENABLED": "OFF",
+                "CMAKE_EXPORT_COMPILE_COMMANDS": "ON",
+                "CMAKE_INSTALL_PREFIX": "${sourceDir}/build/install"
+            },
+            "environment": {
+                "CC": "clang",
+                "CXX": "clang++"
+            },
+            "architecture": {
+                "value": "x64",
+                "strategy": "external"
+            }
+        },
+        {
+            "name": "tests",
+            "inherits": "default",
+            "cacheVariables": {
+                "DSPLIB_BUILD_TESTS": "ON",
+                "DSPLIB_ASAN_ENABLED": "ON"
+            }
+        },
+        {
+            "name": "benchs-float64",
+            "inherits": "default",
+            "cacheVariables": {
+                "DSPLIB_BUILD_BENCHS": "ON",
+                "CMAKE_BUILD_TYPE": "Release"
+            }
+        },
+        {
+            "name": "benchs-float32",
+            "inherits": "default",
+            "cacheVariables": {
+                "DSPLIB_BUILD_BENCHS": "ON",
+                "CMAKE_BUILD_TYPE": "Release",
+                "DSPLIB_USE_FLOAT32": "ON"
+            }
+        },
+        {
+            "name": "examples",
+            "inherits": "default",
+            "cacheVariables": {
+                "DSPLIB_BUILD_EXAMPLES": "ON",
+                "CMAKE_BUILD_TYPE": "Release"
+            }
+        }
+    ]
+}

README.md

@@ -22,7 +22,7 @@ for (int i=0; i < r.size(); ++i) {
     r[i] = x1[i] * x2[i];
 }
 
-auto p = fftw_plan_dft_1d(N, x, spec, FFTW_FORWARD, FFTW_ESTIMATE);
+auto p = fftw_plan_dft_1d(N, x.data(), spec.data(), FFTW_FORWARD, FFTW_ESTIMATE);
 fftw_execute(p);
 fftw_destroy_plan(p);
 ```
@@ -32,7 +32,7 @@ and who likes this:
 ```cpp
 using namespace dsplib;
 x *= 0.3;
-auto power = sum(abs2(*x.slice(lb, rb)));
+auto power = sum(abs2(x.slice(lb, rb)));
 auto r = x1 * x2;
 auto spec = fft(x);
 ```
@@ -75,6 +75,7 @@ arr_cmplx y3 = x1 * x2;
 arr_cmplx y4 = x2 * 1000;
 arr_cmplx y5 = x2.slice(0, 2);
 arr_cmplx y6 = x1 * 2i;
+arr_cmplx y7 = complex(y1); // only explicit conversion
 ```
 
 ### Slicing
@@ -90,18 +91,63 @@ x.slice(-8, 7) ///OUT_OF_RANGE, but numpy returns [0 1 2 3 4 5 6]
 ```
 
 ### Fast Fourier Transform:
-The FFT/IFFT calculation table is cached on first run. To eliminate this behavior, you can use the FftPlan object.
+The `FFT` implementation has no radix size limitations. 
+It supports power-of-two, prime, and semiprime radices. 
+
+If your platform has a faster implementation (e.g., `NE10` on `ARM`), you can set the `DSPLIB_EXCLUDE_FFT=ON` option and implement the `fft_plan_c`, `fft_plan_r`, `ifft_plan_c`, and `ifft_plan_r` functions (see the `FFTW` example).
+
+The tables for the `FFT` are stored in the `LRU` cache and can be recalculated (if the pipeline uses many different bases). Use the `FftPlan` object to avoid this.
 ```cpp
+//FFT fn
 arr_real x = randn(500);
-arr_cmplx y1 = fft(x);  //500
-arr_cmplx y2 = fft(x, 1024); //1024
+arr_cmplx y1 = fft(x);  // real fft, n=500
+arr_cmplx y2 = fft(x, 1024); // real fft, n=1024, zero padding
+arr_cmplx y3 = fft(complex(y1)); // cmplx fft, n=500
+arr_cmplx y4 = rfft(y1); // real fft, equal `fft(x)`
 ```
 
-### Inverse Fast Fourier Transform:
 ```cpp
-arr_cmplx x = 1i * zeros(512);
-x[10] = 1;
-arr_cmplx y = ifft(x);
+//FFT Plan
+const int n = 512;
+const arr_real x = randn(n);
+std::shared_ptr<FftPlanR> plan = fft_plan_r(n);
+
+// real fft, n=512
+arr_cmplx y1 = plan->solve(x); 
+//or
+arr_cmplx y2 = plan->solve(x.slice(0, 512));
+//or
+arr_cmplx y3 = plan->solve(make_span(x.data(), n));
+
+//real fft, n=512, result copy to `r`
+arr_cmplx r(n);
+plan->solve(x, r); 
+//or
+plan->solve(make_span(x.data(), n), make_span(r.data(), n))
+```
+
+
+```cpp
+//IFFT fn
+const int n = 512;
+arr_cmplx x = complex(ones(n));
+arr_cmplx y1 = ifft(x);
+//or
+arr_real y2 = irfft(x.slice(0, n/2+1), n);
+//or
+arr_real y3 = irfft(x);
+```
+
+```cpp
+//IFFT Plan
+const int n = 512;
+const auto x = complex(ones(n));
+auto plan = ifft_plan_r(n);
+
+arr_real y1 = plan->solve(x);
+//or
+arr_real y2;
+plan->solve(make_span(x.data(), n/2+1), make_span(y2.data(), n));
 ```
 
 ### FIR filter:
@@ -217,11 +263,11 @@ auto out = dsplib::resample(in, 2, 1);
 auto out = dsplib::resample(in, 32000, 16000);
 ```
 
-## Building
+## Build
 
 ### Requires:
 - CMake (>=3.10)
-- C++ compiler for C++17 standard (gcc, clang, msvc, mingw)
+- C++17 compiler (exceptions can be disabled)
 
 
 ### Build and install:
@@ -241,7 +287,7 @@ CPMAddPackage(NAME dsplib
     GIT_REPOSITORY 
         "https://github.com/vitalsong/dsplib.git"
     VERSION 
-        0.45.0
+        0.55.3
     OPTIONS
         "DSPLIB_USE_FLOAT32 OFF"
         "DSPLIB_NO_EXCEPTIONS OFF"
@@ -259,6 +305,8 @@ cmake --build build
 ./build/benchs/dsplib-benchs
 ```
 
+### FFT
+
 The implementation of non-power-of-two FFT is based on the general factorization algorithm. It is usually slower, but not critical. 
 
 For prime and semi-prime numbers, the czt algorithm is used, which can be significantly slower (but not as slow as regular DFT).
@@ -310,9 +358,12 @@ BM_KISSFFT/16384/min_time:5.000            98.5 us         98.5 us        69101
 ```
 
 ## TODO:
-- Select FFT backend type (fftw/ne10)
-- Add matrix syntax support
-- Add custom allocator for `base_array<T>` type
-- Add audioread/audiowrite functions (optional libsndfile?)
-- Add chain syntax like `fft(x)->abs2()->pow2db()`
-- Use `const_span<T>` args
+- Add matrix syntax support;
+- Add custom allocator for `base_array<T>` type;
+- Add audioread/audiowrite functions (optional libsndfile?);
+- Add chain syntax like `fft(x)->abs2()->pow2db()`;
+- `SOS` filters;
+- Multichannel resampler;
+- Thread-safe storage for `FFT` (not `thread_local`);
+- Add `chirp`, `conv`, `filter`, `dzt`, `remez` etc.
+- Real/Imag slice for `arr_cmplx`;

benchs/adaptive.cpp

@@ -3,7 +3,7 @@
 
 static std::pair<dsplib::arr_real, dsplib::arr_real> _prepare_frame(int M, int L) {
     const auto h = dsplib::randn(M);
-    auto flt = dsplib::FftFilter(h);
+    auto flt = dsplib::FirFilter(h);
     auto x = dsplib::randn(L);
     auto n = 0.01 * dsplib::randn(L);
     auto d = flt(x) + n;

examples/fftw-backend/fft.cpp

@@ -28,7 +28,7 @@ class FFTWPlanC : public FftPlanC
         fftw_free(in_);
     }
 
-    dsplib::arr_cmplx solve(const dsplib::arr_cmplx& x) const final {
+    dsplib::arr_cmplx solve(span_t<cmplx_t> x) const final {
         DSPLIB_ASSERT(x.size() == n_, "input size must be equal `n`");
         std::memcpy(in_, x.data(), n_ * sizeof(x[0]));
         fftw_execute(plan_);
@@ -66,7 +66,7 @@ class FFTWPlanR : public FftPlanR
         fftw_free(in_);
     }
 
-    dsplib::arr_cmplx solve(const dsplib::arr_real& x) const final {
+    dsplib::arr_cmplx solve(span_t<real_t> x) const final {
         DSPLIB_ASSERT(x.size() == n_, "input size must be equal `n`");
         std::memcpy(in_, x.data(), n_ * sizeof(x[0]));
         fftw_execute(plan_);
@@ -105,7 +105,7 @@ class IFFTWPlanR : public IfftPlanR
         fftw_free(in_);
     }
 
-    dsplib::arr_real solve(const dsplib::arr_cmplx& x) const final {
+    dsplib::arr_real solve(span_t<cmplx_t> x) const final {
         const int n2 = n_ / 2 + 1;
         DSPLIB_ASSERT((x.size() == n_) || (x.size() == n2), "input size must be equal `n` or `n/2+1`");
         std::memcpy(in_, x.data(), n2 * sizeof(x[0]));