C backend for PocketFFT

.github/workflows/ci.yml

@@ -0,0 +1,82 @@
+name: impulse CI
+on:
+  push:
+    paths:
+      - 'tests/**'
+      - '**'
+      - 'impulse.nimble'
+      - '.github/workflows/ci.yml'
+  pull_request:
+    paths:
+      - 'tests/**'
+      - '**'
+      - 'impulse.nimble'
+      - '.github/workflows/ci.yml'
+
+jobs:
+  build:
+    strategy:
+      fail-fast: false
+      matrix:
+        branch: [version-1-6, devel]
+        target: [linux, macos, windows]
+        include:
+          - target: linux
+            builder: ubuntu-latest
+          - target: macos
+            builder: macos-latest
+          - target: windows
+            builder: windows-latest
+    name: '${{ matrix.target }} (${{ matrix.branch }})'
+    runs-on: ${{ matrix.builder }}
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v2
+        with:
+          path: impulse
+
+      - name: Setup Nim
+        uses: alaviss/[email protected]
+        with:
+          path: nim
+          version: ${{ matrix.branch }}
+
+      - name: Setup nimble & deps
+        shell: bash
+        run: |
+          cd impulse
+          nimble refresh -y
+          nimble install -y
+
+      - name: Run tests
+        shell: bash
+        run: |
+          cd impulse
+          nimble -y test
+
+      - name: Build docs
+        if: >
+          github.event_name == 'push' && github.ref == 'refs/heads/master' &&
+          matrix.target == 'linux' && matrix.branch == 'devel'
+        shell: bash
+        run: |
+          cd impulse
+          branch=${{ github.ref }}
+          branch=${branch##*/}
+          nimble doc --project --outdir:docs \
+            '--git.url:https://github.com/${{ github.repository }}' \
+            '--git.commit:${{ github.sha }}' \
+            "--git.devel:$branch" \
+            impulse.nim
+          # Ignore failures for older Nim
+          cp docs/{the,}index.html || true
+
+      - name: Publish docs
+        if: >
+          github.event_name == 'push' && github.ref == 'refs/heads/master' &&
+          matrix.target == 'linux' && matrix.branch == 'devel'
+        uses: crazy-max/ghaction-github-pages@v1
+        with:
+          build_dir: impulse/docs
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

README.md

@@ -5,13 +5,76 @@ Impulse will be a collection of signal processing primitives for Nim.
 ## FFT
 
 Currently this library only consists of an FFT module, which wraps
-[PocketFFT](https://gitlab.mpcdf.mpg.de/mtr/pocketfft) in form of a
-header-only C++ version, https://github.com/mreineck/pocketfft
+[PocketFFT](https://gitlab.mpcdf.mpg.de/mtr/pocketfft).
 
-To use it, please import the submodule directly. For example:
+For the C++ backend an optimized version of PocketFFT is used, in the
+form of a header-only version,
+
+https://github.com/mreineck/pocketfft
+
+For the C backend we use the PocketFFT directly (the C library linked
+before).
+
+Note that the API differs slightly between the two. See the two
+examples below.
+
+### C example
+
+The `pocketfft` submodule can to be imported manually using
+`import impulse/fft/pocketfft` or one can simply import the `fft`
+submodule as shown below.
+
+```nim
+import impulse/fft
+
+template isClose(a1, a2, eps: untyped): untyped =
+  for i in 0 ..< a1.len:
+    doAssert abs(a1[i] - a2[i]) < eps, "Is: " & $a1[i] & " | " & $a2[i]
+
+block Array:
+  let dIn = [1.0, 2.0, 1.0, -1.0, 1.5]
+  let dOut = fft(dIn, forward = true) # forward or backwards?
+  echo dIn
+  echo dOut
+  isClose(dIn, fft(dOut, forward = false), eps = 1e-10)
+  # [1.0, 2.0, 1.0, -1.0, 1.5]
+  # @[4.5, 2.081559480312316, -1.651098762732523, -1.831559480312316, 1.608220406444071]
+block Seq:
+  let dIn = @[1.0, 2.0, 1.0, -1.0, 1.5]
+  let dOut = fft(dIn, forward = true) # forward or backwards?
+  echo dIn
+  echo dOut
+  isClose(dIn, fft(dOut, forward = false), eps = 1e-10)
+  # @[1.0, 2.0, 1.0, -1.0, 1.5]
+  # @[4.5, 2.081559480312316, -1.651098762732523, -1.831559480312316, 1.608220406444071]
+block Tensor:
+  let dIn = @[1.0, 2.0, 1.0, -1.0, 1.5].toTensor
+  let dOut = fft(dIn, forward = true) # forward or backwards?
+  echo dIn
+  echo dOut
+  isClose(dIn, fft(dOut, forward = false), eps = 1e-10)
+  # Tensor[system.float] of shape "[5]" on backend "Cpu"
+  #     1      2      1     -1    1.5
+  # Tensor[system.float] of shape "[5]" on backend "Cpu"
+  #     4.5     2.08156     -1.6511    -1.83156     1.60822
+import std / complex
+block Complex:
+  let dIn = [complex(1.0), complex(2.0), complex(1.0), complex(-1.0), complex(1.5)]
+  let dOut = fft(dIn, forward = true) # forward or backwards?
+  echo dIn
+  echo dOut
+  isClose(dIn, fft(dOut, forward = false), eps = 1e-10)
+  # [(1.0, 0.0), (2.0, 0.0), (1.0, 0.0), (-1.0, 0.0), (1.5, 0.0)]
+  # @[(4.5, 0.0), (2.081559480312316, -1.651098762732523), (-1.831559480312316, 1.608220406444071), (-1.831559480312316, -1.608220406444071), (2.081559480312316, 1.651098762732523)]
+```
+
+
+### C++ example
+
+When compiling on the C++ backend, the API is a bit different:
 
 ```nim
-import impulse/fft/pocketfft
+import impulse/fft
 import std / complex
 
 let dIn = @[1.0, 2.0, 1.0, -1.0, 1.5]

impulse.nimble

@@ -9,3 +9,8 @@ license       = "MIT"
 # Dependencies
 
 requires "nim >= 1.6.0"
+requires "arraymancer >= 0.7.28"
+
+task test, "Run standard tests":
+  exec "nim c -r tests/test_fft.nim"
+  exec "nim c -r tests/test_fft2.nim"

impulse/fft.nim

@@ -4,3 +4,6 @@
 #   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
 #   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
 # at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import fft/pocketfft
+export pocketfft

impulse/fft/README.md

@@ -8,6 +8,10 @@ and in particular uses dynamic memory allocation
 
 We currently use PocketFFT.
 
+For the C++ backend we use a header only version of it. For the C
+backend we use the regular C PocketFFT implementation (which we build
+alongside the Nim program).
+
 Commit 49b813232507470a047727712acda105b84c7815
 has an initial pure Nim implementation that follows
 the PocketFFT algorithms which in turn follow

impulse/fft/c_pocketfft/LICENSE.md

@@ -0,0 +1,25 @@
+Copyright (C) 2010-2019 Max-Planck-Society
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+* Redistributions in binary form must reproduce the above copyright notice, this
+  list of conditions and the following disclaimer in the documentation and/or
+  other materials provided with the distribution.
+* Neither the name of the copyright holder nor the names of its contributors may
+  be used to endorse or promote products derived from this software without
+  specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

impulse/fft/c_pocketfft/README.md

@@ -0,0 +1,55 @@
+PocketFFT
+---------
+
+This is a heavily modified implementation of FFTPack [1,2], with the following
+advantages:
+
+- strictly C99 compliant
+- more accurate twiddle factor computation
+- very fast plan generation
+- worst case complexity for transform sizes with large prime factors is
+  `N*log(N)`, because Bluestein's algorithm [3] is used for these cases.
+
+License
+-------
+
+3-clause BSD (see LICENSE.md)
+
+
+Some code details
+-----------------
+
+Twiddle factor computation:
+
+- making use of symmetries to reduce number of sin/cos evaluations
+- all angles are reduced to the range `[0; pi/4]` for higher accuracy
+- an adapted implementation of `sincospi()` is used, which actually computes
+  `sin(x)` and `(cos(x)-1)`.
+- if `n` sin/cos pairs are required, the adjusted `sincospi()` is only called
+  `2*sqrt(n)` times; the remaining values are obtained by evaluating the
+  angle addition theorems in a numerically accurate way.
+
+Parallel invocation:
+
+- Plans only contain read-only data; all temporary arrays are allocated and
+  deallocated during an individual FFT execution. This means that a single plan
+  can be used in several threads at the same time.
+
+Efficient codelets are available for the factors:
+
+- 2, 3, 4, 5, 7, 11 for complex-valued FFTs
+- 2, 3, 4, 5 for real-valued FFTs
+
+Larger prime factors are handled by somewhat less efficient, generic routines.
+
+For lengths with very large prime factors, Bluestein's algorithm is used, and
+instead of an FFT of length `n`, a convolution of length `n2 >= 2*n-1`
+is performed, where `n2` is chosen to be highly composite.
+
+
+[1] Swarztrauber, P. 1982, Vectorizing the Fast Fourier Transforms
+    (New York: Academic Press), 51
+
+[2] https://www.netlib.org/fftpack/
+
+[3] https://en.wikipedia.org/wiki/Chirp_Z-transform

impulse/fft/c_pocketfft/ffttest.c

@@ -0,0 +1,101 @@
+/*
+ * This file is part of pocketfft.
+ * Licensed under a 3-clause BSD style license - see LICENSE.md
+ */
+
+/*
+ *  Test codes for pocketfft.
+ *
+ *  Copyright (C) 2004-2018 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "pocketfft_src.h"
+
+#define maxlen 10 //8192
+
+static void fill_random (double *data, size_t length)
+  {
+      for (size_t m=0; m<length; ++m){
+	  data[m] = rand()/(RAND_MAX+1.0)-0.5;
+	  printf("Data at i %i = %f\n", m, data[m]);
+      }
+  }
+
+static double errcalc (double *data, double *odata, size_t length)
+  {
+  double sum = 0, errsum = 0;
+  for (size_t m=0; m<length; ++m)
+    {
+    errsum += (data[m]-odata[m])*(data[m]-odata[m]);
+    sum += odata[m]*odata[m];
+    }
+  return sqrt(errsum/sum);
+  }
+
+static int test_real(void)
+  {
+  double data[maxlen], odata[maxlen];
+  const double epsilon=2e-15;
+  int ret = 0;
+  fill_random (odata, maxlen);
+  double errsum=0;
+  for (int length=1; length<=maxlen; ++length)
+    {
+    memcpy (data,odata,length*sizeof(double));
+    rfft_plan plan = make_rfft_plan (length);
+    printf("Input = %f\n", data[0]);
+    rfft_forward (plan, data, 1.);
+    printf("Output = %f\n", data[0]);
+    rfft_backward (plan, data, 1./length);
+    printf("Real output = %f\n", data[0]);
+    destroy_rfft_plan (plan);
+    double err = errcalc (data, odata, length);
+    if (err>epsilon)
+      {
+      printf("problem at real length %i: %e\n",length,err);
+      ret = 1;
+      }
+    errsum+=err;
+    }
+  printf("errsum: %e\n",errsum);
+  return ret;
+  }
+
+static int test_complex(void)
+  {
+  double data[2*maxlen], odata[2*maxlen];
+  fill_random (odata, 2*maxlen);
+  const double epsilon=2e-15;
+  int ret = 0;
+  double errsum=0;
+  for (int length=1; length<=maxlen; ++length)
+    {
+    memcpy (data,odata,2*length*sizeof(double));
+    cfft_plan plan = make_cfft_plan (length);
+    cfft_forward(plan, data, 1.);
+    cfft_backward(plan, data, 1./length);
+    destroy_cfft_plan (plan);
+    double err = errcalc (data, odata, 2*length);
+    if (err>epsilon)
+      {
+      printf("problem at complex length %i: %e\n",length,err);
+      ret = 1;
+      }
+    errsum+=err;
+    }
+  printf("errsum: %e\n",errsum);
+  return ret;
+  }
+
+int main(void)
+  {
+  int ret = 0;
+  ret = test_real();
+  ret += test_complex();
+  return ret;
+  }