Merge pull request #12 from unevens/cpp-wrapper

c++ wrapper with test

3 files changed, 856 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5db02eb..69c794e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -124,6 +124,10 @@ add_executable( test_pffft_double  test_pffft_double.c )
 target_compile_definitions(test_pffft_double PRIVATE _USE_MATH_DEFINES)
 target_link_libraries( test_pffft_double  PFFFT ${ASANLIB} )
 
+add_executable( test_pffft_cpp test_pffft.cpp )
+target_compile_definitions(test_pffft_cpp PRIVATE _USE_MATH_DEFINES)
+target_link_libraries( test_pffft_cpp  PFFFT ${ASANLIB} )
+
 
 add_executable(test_pffastconv   test_pffastconv.c )
 target_compile_definitions(test_pffastconv PRIVATE _USE_MATH_DEFINES)
diff --git a/pffft.hpp b/pffft.hpp
new file mode 100644
index 0000000..75b277e
--- /dev/null
+++ b/pffft.hpp
@@ -0,0 +1,532 @@
+/* Copyright (c) 2013  Julien Pommier ( pommier@modartt.com )
+   Copyright (c) 2020  Hayati Ayguen ( h_ayguen@web.de )
+   Copyright (c) 2020  Dario Mambro ( dario.mambro@gmail.com )
+*/
+
+#pragma once
+
+#include <complex>
+
+namespace {
+#include "pffft.h"
+#include "pffft_double.h"
+}
+
+namespace pffft {
+
+// enum { PFFFT_REAL, PFFFT_COMPLEX }
+typedef pffft_transform_t TransformType;
+
+// enum { PFFFT_FORWARD, PFFFT_BACKWARD }
+typedef pffft_direction_t TransformDirection;
+
+namespace {
+template<typename T>
+class Setup;
+}
+
+// T can be float, double, std::complex<float> or std::complex<double>
+template<typename T>
+class Fft
+{
+public:
+  typedef typename Setup<T>::Scalar Scalar;
+
+  Fft(int length, int stackThresholdLen = 4096);
+
+  ~Fft();
+
+  void prepareLength(int newLength);
+
+  int getLength() const { return length; }
+
+  std::complex<Scalar>* forward(const T* input, std::complex<Scalar>* spectrum);
+
+  T* inverse(const std::complex<Scalar>* spectrum, T* output);
+
+  Scalar* forwardInternalLayout(const T* input,
+                                Scalar* spectrum_internal_layout);
+
+  T* inverseInternalLayout(const Scalar* spectrum_internal_layout, T* output);
+
+  void reorderSpectrum(const Scalar* input,
+                       std::complex<Scalar>* output,
+                       TransformDirection direction);
+
+  Scalar* convolveAccumulate(const Scalar* dft_a,
+                             const Scalar* dft_b,
+                             Scalar* dft_ab,
+                             const Scalar scaling);
+
+  Scalar* convolve(const Scalar* dft_a,
+                   const Scalar* dft_b,
+                   Scalar* dft_ab,
+                   const Scalar scaling);
+
+  template<typename S>
+  static S* alignedAlloc(int length);
+
+  static void alignedFree(void* ptr);
+
+  static Scalar* alignedAllocScalar(int length);
+
+  static std::complex<Scalar>* alignedAllocComplex(int length);
+
+private:
+  Setup<T> setup;
+  Scalar* work;
+  int length;
+  int stackThresholdLen;
+};
+
+// helpers
+
+void
+alignedFree(void* ptr)
+{
+  pffft_aligned_free(ptr);
+}
+
+/* simple helper to get minimum possible fft length */
+int
+minFFtsize(const TransformType transform)
+{
+  return pffft_min_fft_size(transform);
+}
+
+/* simple helper to determine next power of 2
+   - without inexact/rounding floating point operations
+*/
+int
+nextPowerOfTwo(const int N)
+{
+  return pffft_next_power_of_two(N);
+}
+
+/* simple helper to determine if power of 2 - returns bool */
+bool
+isPowerOfTwo(const int N)
+{
+  return pffft_is_power_of_two(N);
+}
+
+// implementation
+
+namespace {
+
+template<typename T>
+struct Setup
+{};
+
+template<>
+class Setup<float>
+{
+  PFFFT_Setup* self;
+
+public:
+  typedef float Scalar;
+
+  Setup()
+    : self(nullptr)
+  {}
+
+  void prepareLength(int length)
+  {
+    if (self) {
+      pffft_destroy_setup(self);
+    }
+    self = pffft_new_setup(length, PFFFT_REAL);
+  }
+
+  ~Setup() { pffft_destroy_setup(self); }
+
+  void transform_ordered(const Scalar* input,
+                         Scalar* output,
+                         Scalar* work,
+                         pffft_direction_t direction)
+  {
+    pffft_transform_ordered(self, input, output, work, direction);
+  }
+
+  void transform(const Scalar* input,
+                 Scalar* output,
+                 Scalar* work,
+                 pffft_direction_t direction)
+  {
+    pffft_transform(self, input, output, work, direction);
+  }
+
+  void reorder(const Scalar* input, Scalar* output, pffft_direction_t direction)
+  {
+    pffft_zreorder(self, input, output, direction);
+  }
+
+  void convolveAccumulate(const Scalar* dft_a,
+                          const Scalar* dft_b,
+                          Scalar* dft_ab,
+                          const Scalar scaling)
+  {
+    pffft_zconvolve_accumulate(self, dft_a, dft_b, dft_ab, scaling);
+  }
+
+  void convolve(const Scalar* dft_a,
+                const Scalar* dft_b,
+                Scalar* dft_ab,
+                const Scalar scaling)
+  {
+    pffft_zconvolve_no_accu(self, dft_a, dft_b, dft_ab, scaling);
+  }
+
+  template<typename S>
+  static S* allocate(int length)
+  {
+    const int bytes = sizeof(S) * length;
+    return static_cast<S*>(pffft_aligned_malloc(bytes));
+  }
+};
+
+template<>
+class Setup<std::complex<float>>
+{
+  PFFFT_Setup* self;
+
+public:
+  typedef float Scalar;
+
+  Setup()
+    : self(nullptr)
+  {}
+
+  ~Setup() { pffft_destroy_setup(self); }
+
+  void prepareLength(int length)
+  {
+    if (self) {
+      pffft_destroy_setup(self);
+    }
+    self = pffft_new_setup(length, PFFFT_COMPLEX);
+  }
+
+  void transform_ordered(const Scalar* input,
+                         Scalar* output,
+                         Scalar* work,
+                         pffft_direction_t direction)
+  {
+    pffft_transform_ordered(self, input, output, work, direction);
+  }
+
+  void transform(const Scalar* input,
+                 Scalar* output,
+                 Scalar* work,
+                 pffft_direction_t direction)
+  {
+    pffft_transform(self, input, output, work, direction);
+  }
+
+  void reorder(const Scalar* input, Scalar* output, pffft_direction_t direction)
+  {
+    pffft_zreorder(self, input, output, direction);
+  }
+
+  void convolve(const Scalar* dft_a,
+                const Scalar* dft_b,
+                Scalar* dft_ab,
+                const Scalar scaling)
+  {
+    pffft_zconvolve_no_accu(self, dft_a, dft_b, dft_ab, scaling);
+  }
+
+  template<typename S>
+  static S* allocate(const int length)
+  {
+    const int bytes = sizeof(S) * length;
+    return static_cast<S*>(pffft_aligned_malloc(bytes));
+  }
+};
+
+template<>
+class Setup<double>
+{
+  PFFFTD_Setup* self;
+
+public:
+  typedef double Scalar;
+
+  Setup()
+    : self(nullptr)
+  {}
+
+  ~Setup() { pffftd_destroy_setup(self); }
+
+  void prepareLength(int length)
+  {
+    if (self) {
+      pffftd_destroy_setup(self);
+      self = nullptr;
+    }
+    if (length > 0) {
+      self = pffftd_new_setup(length, PFFFT_REAL);
+    }
+  }
+
+  void transform_ordered(const Scalar* input,
+                         Scalar* output,
+                         Scalar* work,
+                         pffft_direction_t direction)
+  {
+    pffftd_transform_ordered(self, input, output, work, direction);
+  }
+
+  void transform(const Scalar* input,
+                 Scalar* output,
+                 Scalar* work,
+                 pffft_direction_t direction)
+  {
+    pffftd_transform(self, input, output, work, direction);
+  }
+
+  void reorder(const Scalar* input, Scalar* output, pffft_direction_t direction)
+  {
+    pffftd_zreorder(self, input, output, direction);
+  }
+
+  void convolveAccumulate(const Scalar* dft_a,
+                          const Scalar* dft_b,
+                          Scalar* dft_ab,
+                          const Scalar scaling)
+  {
+    pffftd_zconvolve_accumulate(self, dft_a, dft_b, dft_ab, scaling);
+  }
+
+  void convolve(const Scalar* dft_a,
+                const Scalar* dft_b,
+                Scalar* dft_ab,
+                const Scalar scaling)
+  {
+    pffftd_zconvolve_no_accu(self, dft_a, dft_b, dft_ab, scaling);
+  }
+
+  template<typename S>
+  static S* allocate(int length)
+  {
+    const int bytes = sizeof(S) * length;
+    return static_cast<S*>(pffftd_aligned_malloc(bytes));
+  }
+};
+
+template<>
+class Setup<std::complex<double>>
+{
+  PFFFTD_Setup* self;
+
+public:
+  typedef double Scalar;
+
+  Setup()
+    : self(nullptr)
+  {}
+
+  ~Setup() { pffftd_destroy_setup(self); }
+
+  void prepareLength(int length)
+  {
+    if (self) {
+      pffftd_destroy_setup(self);
+    }
+    self = pffftd_new_setup(length, PFFFT_COMPLEX);
+  }
+
+  void transform_ordered(const Scalar* input,
+                         Scalar* output,
+                         Scalar* work,
+                         pffft_direction_t direction)
+  {
+    pffftd_transform_ordered(self, input, output, work, direction);
+  }
+
+  void transform(const Scalar* input,
+                 Scalar* output,
+                 Scalar* work,
+                 pffft_direction_t direction)
+  {
+    pffftd_transform(self, input, output, work, direction);
+  }
+
+  void reorder(const Scalar* input, Scalar* output, pffft_direction_t direction)
+  {
+    pffftd_zreorder(self, input, output, direction);
+  }
+
+  void convolveAccumulate(const Scalar* dft_a,
+                          const Scalar* dft_b,
+                          Scalar* dft_ab,
+                          const Scalar scaling)
+  {
+    pffftd_zconvolve_accumulate(self, dft_a, dft_b, dft_ab, scaling);
+  }
+
+  void convolve(const Scalar* dft_a,
+                const Scalar* dft_b,
+                Scalar* dft_ab,
+                const Scalar scaling)
+  {
+    pffftd_zconvolve_no_accu(self, dft_a, dft_b, dft_ab, scaling);
+  }
+
+  template<typename S>
+  static S* allocate(int length)
+  {
+    const int bytes = sizeof(S) * length;
+    return static_cast<S*>(pffftd_aligned_malloc(bytes));
+  }
+};
+
+}
+
+template<typename T>
+inline Fft<T>::Fft(int length, int stackThresholdLen)
+  : length(0)
+  , stackThresholdLen(stackThresholdLen)
+  , work(nullptr)
+{
+  prepareLength(length);
+}
+
+template<typename T>
+inline Fft<T>::~Fft()
+{
+  pffft_aligned_free(work);
+}
+
+template<typename T>
+inline void
+Fft<T>::prepareLength(int newLength)
+{
+  const bool wasOnHeap = work != nullptr;
+
+  const bool useHeap = newLength > stackThresholdLen;
+
+  if (useHeap == wasOnHeap && newLength == length) {
+    return;
+  }
+
+  length = newLength;
+
+  setup.prepareLength(length);
+
+  if (work) {
+    pffft_aligned_free(work);
+    work = nullptr;
+  }
+
+  if (useHeap) {
+    int const bytesToAllocate = length * sizeof(T);
+    work = static_cast<Scalar*>(pffft_aligned_malloc(bytesToAllocate));
+  }
+}
+
+template<typename T>
+inline std::complex<typename Fft<T>::Scalar>*
+Fft<T>::forward(const T* input, std::complex<Scalar>* spectrum)
+{
+  setup.transform_ordered(reinterpret_cast<const Scalar*>(input),
+                          reinterpret_cast<Scalar*>(spectrum),
+                          work,
+                          PFFFT_FORWARD);
+  return spectrum;
+}
+
+template<typename T>
+inline T*
+Fft<T>::inverse(std::complex<Scalar> const* spectrum, T* output)
+{
+  setup.transform_ordered(reinterpret_cast<const Scalar*>(spectrum),
+                          reinterpret_cast<Scalar*>(output),
+                          work,
+                          PFFFT_BACKWARD);
+  return output;
+}
+
+template<typename T>
+inline typename pffft::Fft<T>::Scalar*
+Fft<T>::forwardInternalLayout(const T* input, Scalar* spectrum_internal_layout)
+{
+  setup.transform(reinterpret_cast<const Scalar*>(input),
+                  spectrum_internal_layout,
+                  work,
+                  PFFFT_FORWARD);
+  return spectrum_internal_layout;
+}
+
+template<typename T>
+inline T*
+Fft<T>::inverseInternalLayout(const Scalar* spectrum_internal_layout, T* output)
+{
+  setup.transform(spectrum_internal_layout,
+                  reinterpret_cast<Scalar*>(output),
+                  work,
+                  PFFFT_BACKWARD);
+  return output;
+}
+
+template<typename T>
+inline void
+Fft<T>::reorderSpectrum(const Scalar* input,
+                        std::complex<Scalar>* output,
+                        TransformDirection direction)
+{
+  setup.reorder(input, reinterpret_cast<Scalar*>(output), direction);
+}
+
+template<typename T>
+inline typename pffft::Fft<T>::Scalar*
+Fft<T>::convolveAccumulate(const Scalar* dft_a,
+                           const Scalar* dft_b,
+                           Scalar* dft_ab,
+                           const Scalar scaling)
+{
+  setup.convolveAccumulate(dft_a, dft_b, dft_ab, scaling);
+  return dft_ab;
+}
+
+template<typename T>
+inline typename pffft::Fft<T>::Scalar*
+Fft<T>::convolve(const Scalar* dft_a,
+                 const Scalar* dft_b,
+                 Scalar* dft_ab,
+                 const Scalar scaling)
+{
+  setup.convolve(dft_a, dft_b, dft_ab, scaling);
+  return dft_ab;
+}
+
+template<typename T>
+inline void
+Fft<T>::alignedFree(void* ptr)
+{
+  pffft_aligned_free(ptr);
+}
+
+template<typename T>
+inline typename pffft::Fft<T>::Scalar*
+pffft::Fft<T>::alignedAllocScalar(int length)
+{
+  return alignedAlloc<Scalar>(length);
+}
+
+template<typename T>
+inline std::complex<typename pffft::Fft<T>::Scalar>*
+Fft<T>::alignedAllocComplex(int length)
+{
+  return alignedAlloc<std::complex<Scalar>>(length);
+}
+
+template<typename T>
+template<typename S>
+inline S*
+Fft<T>::alignedAlloc(int length)
+{
+  return Setup<T>::allocate<S>(length);
+}
+
+} // namespace pffft
diff --git a/test_pffft.cpp b/test_pffft.cpp
new file mode 100644
index 0000000..050a488
--- /dev/null
+++ b/test_pffft.cpp
@@ -0,0 +1,320 @@
+/*
+  Copyright (c) 2013 Julien Pommier.
+
+  Small test & bench for PFFFT, comparing its performance with the scalar
+  FFTPACK, FFTW, and Apple vDSP
+
+  How to build:
+
+  on linux, with fftw3:
+  gcc -o test_pffft -DHAVE_FFTW -msse -mfpmath=sse -O3 -Wall -W pffft.c
+  test_pffft.c fftpack.c -L/usr/local/lib -I/usr/local/include/ -lfftw3f -lm
+
+  on macos, without fftw3:
+  clang -o test_pffft -DHAVE_VECLIB -O3 -Wall -W pffft.c test_pffft.c fftpack.c
+  -L/usr/local/lib -I/usr/local/include/ -framework Accelerate
+
+  on macos, with fftw3:
+  clang -o test_pffft -DHAVE_FFTW -DHAVE_VECLIB -O3 -Wall -W pffft.c
+  test_pffft.c fftpack.c -L/usr/local/lib -I/usr/local/include/ -lfftw3f
+  -framework Accelerate
+
+  as alternative: replace clang by gcc.
+
+  on windows, with visual c++:
+  cl /Ox -D_USE_MATH_DEFINES /arch:SSE test_pffft.c pffft.c fftpack.c
+
+  build without SIMD instructions:
+  gcc -o test_pffft -DPFFFT_SIMD_DISABLE -O3 -Wall -W pffft.c test_pffft.c
+  fftpack.c -lm
+
+ */
+
+#include "pffft.hpp"
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+
+/* maximum allowed phase error in degree */
+#define DEG_ERR_LIMIT 1E-4
+
+/* maximum allowed magnitude error in amplitude (of 1.0 or 1.1) */
+#define MAG_ERR_LIMIT 1E-6
+
+#define PRINT_SPEC 0
+
+#define PWR2LOG(PWR) ((PWR) < 1E-30 ? 10.0 * log10(1E-30) : 10.0 * log10(PWR))
+
+template<typename T>
+bool
+Ttest(int N, bool useOrdered)
+{
+  using Fft = typename pffft::Fft<T>;
+  using Scalar = typename Fft::Scalar;
+
+  bool cplx = std::is_same<T, std::complex<float>>::value ||
+              std::is_same<T, std::complex<double>>::value;
+
+  double EXPECTED_DYN_RANGE =
+    std::is_same<double, Scalar>::value ? 215.0 : 140.0;
+
+  int Nsca = (cplx ? N * 2 : N);
+  int Ncplx = (cplx ? N : N / 2);
+  T* X = Fft::alignedAlloc<T>(Nsca);
+  T* Z = Fft::alignedAlloc<T>(Nsca);
+  Scalar* R = Fft::alignedAllocScalar(Nsca);
+  std::complex<Scalar>* Y = Fft::alignedAllocComplex(Nsca);
+  int k, j, m, iter, kmaxOther;
+  bool retError = false;
+  double freq, dPhi, phi, phi0;
+  double pwr, pwrCar, pwrOther, err, errSum, mag, expextedMag;
+  double amp = 1.0;
+
+  assert(pffft::isPowerOfTwo(N));
+
+  Fft fft = Fft(N);
+
+  Scalar* Xs = reinterpret_cast<Scalar*>(X);
+  Scalar* Ys = reinterpret_cast<Scalar*>(Y);
+  Scalar* Zs = reinterpret_cast<Scalar*>(Z);
+
+  for (k = m = 0; k < (cplx ? N : (1 + N / 2)); k += N / 16, ++m) {
+    amp = ((m % 3) == 0) ? 1.0F : 1.1F;
+    freq = (k < N / 2) ? ((double)k / N) : ((double)(k - N) / N);
+    dPhi = 2.0 * M_PI * freq;
+    if (dPhi < 0.0)
+      dPhi += 2.0 * M_PI;
+
+    iter = -1;
+    while (1) {
+      ++iter;
+
+      if (iter)
+        printf("bin %d: dphi = %f for freq %f\n", k, dPhi, freq);
+
+      /* generate cosine carrier as time signal - start at defined phase phi0 */
+      phi = phi0 =
+        (m % 4) * 0.125 * M_PI; /* have phi0 < 90 deg to be normalized */
+      for (j = 0; j < N; ++j) {
+        if (cplx) {
+          Xs[2 * j] = amp * cos(phi);     /* real part */
+          Xs[2 * j + 1] = amp * sin(phi); /* imag part */
+        } else
+          Xs[j] = amp * cos(phi); /* only real part */
+
+        /* phase increment .. stay normalized - cos()/sin() might degrade! */
+        phi += dPhi;
+        if (phi >= M_PI)
+          phi -= 2.0 * M_PI;
+      }
+
+      /* forward transform from X --> Y  .. using work buffer W */
+      if (useOrdered)
+        fft.forward(X, Y);
+      else {
+        fft.forwardInternalLayout(X, R); /* temporarily use R for reordering */
+        fft.reorderSpectrum(R, Y, PFFFT_FORWARD);
+      }
+
+      pwrOther = -1.0;
+      pwrCar = 0;
+
+      /* for positive frequencies: 0 to 0.5 * samplerate */
+      /* and also for negative frequencies: -0.5 * samplerate to 0 */
+      for (j = 0; j < (cplx ? N : (1 + N / 2)); ++j) {
+        if (!cplx && !j) /* special treatment for DC for real input */
+          pwr = Ys[j] * Ys[j];
+        else if (!cplx && j == N / 2) /* treat 0.5 * samplerate */
+          pwr = Ys[1] *
+                Ys[1]; /* despite j (for freq calculation) we have index 1 */
+        else
+          pwr = Ys[2 * j] * Ys[2 * j] + Ys[2 * j + 1] * Ys[2 * j + 1];
+        if (iter || PRINT_SPEC)
+          printf("%s fft %d:  pwr[j = %d] = %g == %f dB\n",
+                 (cplx ? "cplx" : "real"),
+                 N,
+                 j,
+                 pwr,
+                 PWR2LOG(pwr));
+        if (k == j)
+          pwrCar = pwr;
+        else if (pwr > pwrOther) {
+          pwrOther = pwr;
+          kmaxOther = j;
+        }
+      }
+
+      if (PWR2LOG(pwrCar) - PWR2LOG(pwrOther) < EXPECTED_DYN_RANGE) {
+        printf("%s fft %d amp %f iter %d:\n",
+               (cplx ? "cplx" : "real"),
+               N,
+               amp,
+               iter);
+        printf("  carrier power  at bin %d: %g == %f dB\n",
+               k,
+               pwrCar,
+               PWR2LOG(pwrCar));
+        printf("  carrier mag || at bin %d: %g\n", k, sqrt(pwrCar));
+        printf("  max other pwr  at bin %d: %g == %f dB\n",
+               kmaxOther,
+               pwrOther,
+               PWR2LOG(pwrOther));
+        printf("  dynamic range: %f dB\n\n",
+               PWR2LOG(pwrCar) - PWR2LOG(pwrOther));
+        retError = true;
+        if (iter == 0)
+          continue;
+      }
+
+      if (k > 0 && k != N / 2) {
+        phi = atan2(Ys[2 * k + 1], Ys[2 * k]);
+        if (fabs(phi - phi0) > DEG_ERR_LIMIT * M_PI / 180.0) {
+          retError = true;
+          printf("%s fft %d  bin %d amp %f : phase mismatch! phase = %f deg   "
+                 "expected = %f deg\n",
+                 (cplx ? "cplx" : "real"),
+                 N,
+                 k,
+                 amp,
+                 phi * 180.0 / M_PI,
+                 phi0 * 180.0 / M_PI);
+        }
+      }
+
+      expextedMag = cplx ? amp : ((k == 0 || k == N / 2) ? amp : (amp / 2));
+      mag = sqrt(pwrCar) / N;
+      if (fabs(mag - expextedMag) > MAG_ERR_LIMIT) {
+        retError = true;
+        printf("%s fft %d  bin %d amp %f : mag = %g   expected = %g\n",
+               (cplx ? "cplx" : "real"),
+               N,
+               k,
+               amp,
+               mag,
+               expextedMag);
+      }
+
+      /* now convert spectrum back */
+      fft.inverse(Y, Z);
+
+      errSum = 0.0;
+      for (j = 0; j < (cplx ? (2 * N) : N); ++j) {
+        /* scale back */
+        Z[j] /= N;
+        /* square sum errors over real (and imag parts) */
+        err = (Xs[j] - Zs[j]) * (Xs[j] - Zs[j]);
+        errSum += err;
+      }
+
+      if (errSum > N * 1E-7) {
+        retError = true;
+        printf("%s fft %d  bin %d : inverse FFT doesn't match original signal! "
+               "errSum = %g ; mean err = %g\n",
+               (cplx ? "cplx" : "real"),
+               N,
+               k,
+               errSum,
+               errSum / N);
+      }
+
+      break;
+    }
+  }
+  pffft::alignedFree(X);
+  pffft::alignedFree(Y);
+  pffft::alignedFree(Z);
+
+  return retError;
+}
+
+bool
+test(int N, bool useComplex, bool useOrdered)
+{
+  if (useComplex) {
+    return Ttest<std::complex<float>>(N, useOrdered) &&
+           Ttest<std::complex<double>>(N, useOrdered);
+  } else {
+    return Ttest<float>(N, useOrdered) && Ttest<double>(N, useOrdered);
+  }
+}
+
+int
+main(int argc, char** argv)
+{
+  int N, result, resN, resAll, k, resNextPw2, resIsPw2, resFFT;
+
+  int inp_power_of_two[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 511, 512, 513 };
+  int ref_power_of_two[] = { 1, 2, 4, 4, 8, 8, 8, 8, 16, 512, 512, 1024 };
+
+  resNextPw2 = 0;
+  resIsPw2 = 0;
+  for (k = 0; k < (sizeof(inp_power_of_two) / sizeof(inp_power_of_two[0]));
+       ++k) {
+    N = pffft::nextPowerOfTwo(inp_power_of_two[k]);
+    if (N != ref_power_of_two[k]) {
+      resNextPw2 = 1;
+      printf("pffft_next_power_of_two(%d) does deliver %d, which is not "
+             "reference result %d!\n",
+             inp_power_of_two[k],
+             N,
+             ref_power_of_two[k]);
+    }
+
+    result = pffft::isPowerOfTwo(inp_power_of_two[k]);
+    if (inp_power_of_two[k] == ref_power_of_two[k]) {
+      if (!result) {
+        resIsPw2 = 1;
+        printf("pffft_is_power_of_two(%d) delivers false; expected true!\n",
+               inp_power_of_two[k]);
+      }
+    } else {
+      if (result) {
+        resIsPw2 = 1;
+        printf("pffft_is_power_of_two(%d) delivers true; expected false!\n",
+               inp_power_of_two[k]);
+      }
+    }
+  }
+  if (!resNextPw2)
+    printf("tests for pffft_next_power_of_two() succeeded successfully.\n");
+  if (!resIsPw2)
+    printf("tests for pffft_is_power_of_two() succeeded successfully.\n");
+
+  resFFT = 0;
+  for (N = 32; N <= 65536; N *= 2) {
+    result = test(N, 1 /* cplx fft */, 1 /* useOrdered */);
+    resN = result;
+    resFFT |= result;
+
+    result = test(N, 0 /* cplx fft */, 1 /* useOrdered */);
+    resN |= result;
+    resFFT |= result;
+
+    result = test(N, 1 /* cplx fft */, 0 /* useOrdered */);
+    resN |= result;
+    resFFT |= result;
+
+    result = test(N, 0 /* cplx fft */, 0 /* useOrdered */);
+    resN |= result;
+    resFFT |= result;
+
+    if (!resN)
+      printf("tests for size %d succeeded successfully.\n", N);
+  }
+
+  if (!resFFT)
+    printf("all pffft transform tests (FORWARD/BACKWARD, REAL/COMPLEX) "
+           "succeeded successfully.\n");
+
+  resAll = resNextPw2 | resIsPw2 | resFFT;
+  if (!resAll)
+    printf("all tests succeeded successfully.\n");
+  else
+    printf("there are failed tests!\n");
+
+  return resAll;
+}