First work on avpowerspectra implementation. Added stubs for all other code that needs to be implemented in C++ as well.

.gitmodules

@@ -16,3 +16,6 @@
 [submodule "third_party/lockfree"]
 	path = third_party/lockfree
 	url = https://github.com/boostorg/lockfree
+[submodule "third_party/carma"]
+	path = third_party/carma
+	url = https://github.com/RUrlus/carma

CMakeLists.txt

@@ -26,7 +26,7 @@ option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ON)
 option(LASP_HAS_ULDAQ "Compile with UlDaq backend" ON)
 
 set(LASP_FFT_BACKEND "FFTW" CACHE STRING "FFT Library backend")
-set_property(CACHE LASP_FFT_BACKEND PROPERTY STRINGS "FFTW" "FFTPack" "None")
+set_property(CACHE LASP_FFT_BACKEND PROPERTY STRINGS "FFTW" "Armadillo")
 
 set(PY_FULL_VERSION ${PROJECT_VERSION}${PY_VERSION_SUFFIX})
 
@@ -52,10 +52,16 @@ include(OSSpecific)
 
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-type-limits")
 
-set(CMAKE_C_FLAGS_RELEASE  "-O3 -mfpmath=sse -march=x86-64 -mtune=native \
+set(CMAKE_C_FLAGS_RELEASE  "-O3 -flto -mfpmath=sse -march=x86-64 -mtune=native \
 -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions")
 
 # ############################# End compilation flags
+add_subdirectory(third_party/carma)
 
-add_subdirectory(third_party)
+if(LASP_FFT_BACKEND STREQUAL "FFTW")
+    find_library(FFTW_LIBRARY NAMES fftw3 fftw)
+    set(FFT_LIBRARIES "${FFTW_LIBRARY}")
+elseif(LASP_FFT_BACKEND STREQUAL "Armadillo")
+endif()
 add_subdirectory(src/lasp)
+

src/lasp/__init__.py

@@ -16,7 +16,7 @@ from .lasp_measurement import *  # Measurement, scaleBlockSens
 from .lasp_octavefilter import *
 # from .lasp_slm import *  # SLM, Dummy
 from .lasp_record import *  # RecordStatus, Recording
-from .lasp_daqconfigs import DaqConfigurations
+from .lasp_daqconfigs import *
 # from .lasp_siggen import *  # SignalType, NoiseType, SiggenMessage, SiggenData, Siggen
 # from .lasp_weighcal import *  # WeighCal
 # from .tools import *  # SmoothingType, smoothSpectralData, SmoothingWidth

src/lasp/c/lasp_fft.c

@@ -1,263 +0,0 @@
-// lasp_fft.c
-//
-// Author: J.A. de Jong - ASCEE
-// 
-// Description:
-// 
-//////////////////////////////////////////////////////////////////////
-#define TRACERPLUS (-5)
-#include "lasp_types.h"
-#include "lasp_config.h"
-#include "lasp_tracer.h"
-#include "lasp_fft.h"
-
-#if LASP_FFT_BACKEND == FFTPack
-#include "fftpack.h"
-typedef struct Fft_s {
-    us nfft;
-    vd fft_work; // Storage memory for fftpack
-} Fft_s;
-#elif LASP_FFT_BACKEND == FFTW
-#include <fftw3.h>
-
-typedef struct Fft_s {
-    us nfft;
-    fftw_plan forward_plan;
-    fftw_plan reverse_plan;
-    c* complex_storage;
-    d* real_storage;
-} Fft_s;
-#else
-#error "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW"
-#endif
-
-void load_fft_wisdom(const char* wisdom) {
-#if LASP_FFT_BACKEND == FFTPack
-#elif LASP_FFT_BACKEND == FFTW
-    if(wisdom) {
-        int rv= fftw_import_wisdom_from_string(wisdom);
-        if(rv != 1) {
-            fprintf(stderr, "Error loading FFTW wisdom");
-        }
-    }
-#endif
-}
-
-char* store_fft_wisdom() {
-#if LASP_FFT_BACKEND == FFTPack
-    return NULL;
-#elif LASP_FFT_BACKEND == FFTW
-    return fftw_export_wisdom_to_string();
-#endif
-}
-
-Fft* Fft_create(const us nfft) {
-    fsTRACE(15);
-    if(nfft == 0) {
-        WARN("nfft should be > 0");
-        return NULL;
-    }
-
-    Fft* fft = a_malloc(sizeof(Fft));
-
-    fft->nfft = nfft;
-
-#if LASP_FFT_BACKEND == FFTPack
-    /* Initialize foreign fft lib */
-    fft->fft_work = vd_alloc(2*nfft+15);
-    npy_rffti(nfft,getvdval(&fft->fft_work,0));
-    check_overflow_vx(fft->fft_work);
-#elif LASP_FFT_BACKEND == FFTW
-    fft->complex_storage = fftw_malloc(sizeof(c) * (nfft/2 + 1));
-    fft->real_storage = fftw_malloc(sizeof(d) * nfft);
-
-    fft->forward_plan = fftw_plan_dft_r2c_1d(nfft,
-            fft->real_storage,
-            fft->complex_storage,
-            FFTW_MEASURE);
-    fft->reverse_plan = fftw_plan_dft_c2r_1d(nfft,
-            fft->complex_storage,
-            fft->real_storage,
-            FFTW_MEASURE);
-
-#endif
-
-    /* print_vd(&fft->fft_work); */
-
-    feTRACE(15); 
-    return fft;
-}
-void Fft_free(Fft* fft) {
-    fsTRACE(15);
-    dbgassert(fft,NULLPTRDEREF);
-#if LASP_FFT_BACKEND == FFTPack
-    vd_free(&fft->fft_work);
-#elif LASP_FFT_BACKEND == FFTW
-    fftw_free(fft->complex_storage);
-    fftw_free(fft->real_storage);
-    fftw_destroy_plan(fft->forward_plan);
-    fftw_destroy_plan(fft->reverse_plan);
-#endif
-    a_free(fft);
-    feTRACE(15);
-}
-
-us Fft_nfft(const Fft* fft) {return fft->nfft;}
-
-void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) {
-    fsTRACE(15);
-    dbgassert(fft && freqdata && result,NULLPTRDEREF);
-    const us nfft = fft->nfft;
-    dbgassert(result->n_rows == nfft,
-            "Invalid size for time data rows."
-            " Should be equal to nfft");
-
-    dbgassert(freqdata->n_rows == (nfft/2+1),"Invalid number of rows in"
-            " result array");
-
-
-    d* result_ptr = getvdval(result,0);
-
-#if LASP_FFT_BACKEND == FFTPack
-    d* freqdata_ptr = (d*) getvcval(freqdata,0);
-    /* Copy freqdata, to fft_result. */
-    d_copy(&result_ptr[1],&freqdata_ptr[2],nfft-1,1,1);
-    result_ptr[0] = freqdata_ptr[0];
-
-    /* Perform inplace backward transform */
-    npy_rfftb(nfft,
-            result_ptr,
-            getvdval(&fft->fft_work,0));
-
-
-#elif LASP_FFT_BACKEND == FFTW
-    c* freqdata_ptr = (c*) getvcval(freqdata,0);
-
-    c_copy(fft->complex_storage, freqdata_ptr,nfft/2+1);
-
-    fftw_execute(fft->reverse_plan);
-    
-    d_copy(result_ptr, fft->real_storage, nfft, 1, 1);
-    
-#endif
-    check_overflow_vx(*result);
-
-    /* Scale by dividing by nfft. Checked with numpy implementation
-     * that this indeed needs to be done for FFTpack. */
-    d_scale(result_ptr,1/((d) nfft),nfft);
-    feTRACE(15);
-}
-void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata) {
-    fsTRACE(15);
-
-    dbgassert(fft && timedata && freqdata,NULLPTRDEREF);
-
-    const us nchannels = timedata->n_cols;
-    dbgassert(timedata->n_cols == freqdata->n_cols,
-            "Number of columns in timedata and result"
-            " should be equal.");
-
-    for(us col=0;col<nchannels;col++) {
-
-        vd timedata_col = dmat_column(timedata,col);
-        vc freqdata_col = cmat_column((cmat*)freqdata,col);
-
-        Fft_ifft_single(fft,&freqdata_col,&timedata_col);
-
-        vd_free(&timedata_col);
-        vc_free(&freqdata_col);
-    }
-    check_overflow_xmat(*timedata);
-    check_overflow_xmat(*freqdata);    
-
-    feTRACE(15);
-}
-
-void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
-
-    fsTRACE(15);
-    dbgassert(fft && timedata && result,NULLPTRDEREF);
-
-    const us nfft = fft->nfft;
-    assert_vx(timedata);
-    assert_vx(result);
-    dbgassert(timedata->n_rows == nfft,
-            "Invalid size for time data rows."
-            " Should be equal to nfft");
-
-    dbgassert(result->n_rows == (nfft/2+1),"Invalid number of rows in"
-            " result array");
-
-
-#if LASP_FFT_BACKEND == FFTPack
-    d* result_ptr = (d*) getvcval(result,0);
-
-    /* Fftpack stores the data a bit strange, the resulting array
-     * has the DC value at 0,the first cosine at 1, the first sine
-     * at 2 etc. 1
-     * resulting matrix, as for the complex data, the imaginary
-     * part of the DC component equals zero. */
-
-    /* Copy timedata, as it will be overwritten in the fft pass. */
-    d_copy(&result_ptr[1],getvdval(timedata,0),nfft,1,1);
-
-
-    /* Perform fft */
-    npy_rfftf(nfft,&result_ptr[1],
-            getvdval(&fft->fft_work,0));
-
-    /* Set real part of DC component to first index of the rfft
-     * routine */
-    result_ptr[0] = result_ptr[1];
-
-    result_ptr[1] = 0;        /* Set imaginary part of DC component
-                               * to zero */
-
-    /* For an even fft, the imaginary part of the Nyquist frequency
-     * bin equals zero.*/
-    if(islikely(nfft%2 == 0)) {
-        result_ptr[nfft+1] = 0;
-    }
-    check_overflow_vx(fft->fft_work);
-#elif LASP_FFT_BACKEND == FFTW
-
-    d* timedata_ptr = getvdval(timedata,0);
-    c* result_ptr = getvcval(result,0);
-    d_copy(fft->real_storage,timedata_ptr, nfft, 1, 1);
-
-    fftw_execute(fft->forward_plan);
-    
-    c_copy(result_ptr, fft->complex_storage, nfft/2+1);
-    
-#endif
-    check_overflow_vx(*result);
-    feTRACE(15);
-
-}
-void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result) {
-    fsTRACE(15);
-
-    dbgassert(fft && timedata && result,NULLPTRDEREF);
-
-    const us nchannels = timedata->n_cols;
-    dbgassert(timedata->n_cols == result->n_cols,
-            "Number of columns in timedata and result"
-            " should be equal.");
-
-    for(us col=0;col<nchannels;col++) {
-
-        vd timedata_col = dmat_column((dmat*) timedata,col);
-        vc result_col = cmat_column(result,col);
-
-        Fft_fft_single(fft,&timedata_col,&result_col);
-
-        vd_free(&timedata_col);
-        vc_free(&result_col);
-    }
-    check_overflow_xmat(*timedata);
-    check_overflow_xmat(*result);
-
-    feTRACE(15);
-}
-
-//////////////////////////////////////////////////////////////////////

src/lasp/device/lasp_daq.h

@@ -4,7 +4,6 @@
 #include "lasp_deviceinfo.h"
 #include "lasp_types.h"
 #include <functional>
-#include <gsl/gsl-lite.hpp>
 #include <memory>
 
 /**

src/lasp/dsp/CMakeLists.txt

@@ -5,6 +5,8 @@ set(lasp_dsp_files
     lasp_siggen.cpp
     lasp_siggen_impl.cpp
     lasp_window.cpp
+    lasp_fft.cpp
+    lasp_avpowerspectra.cpp
     )
 
 

src/lasp/dsp/lasp_avpowerspectra.cpp

@@ -0,0 +1,68 @@
+#define DEBUGTRACE_ENABLED
+#include "lasp_avpowerspectra.h"
+#include "debugtrace.hpp"
+#include <cmath>
+
+using std::runtime_error;
+
+PowerSpectra::PowerSpectra(const us nfft, const Window::WindowType w)
+  : PowerSpectra(Window::create(w, nfft)) {}
+
+PowerSpectra::PowerSpectra(const vd &window)
+  : nfft(window.size()), _fft(nfft), _window(window) {
+
+    d win_pow = arma::sum(window % window);
+
+    /* Scale fft such that power is easily computed */
+    _scale_fac = sqrt(2 / win_pow) / nfft;
+  }
+arma::Cube<c> PowerSpectra::compute(const dmat &input) {
+
+  dmat input_tmp = input;
+
+  // Multiply each column of the inputs element-wise with the window.
+  input_tmp.each_col() %= _window;
+
+  cmat rfft = _fft.fft(input_tmp) * _scale_fac;
+
+  arma::cx_cube output(rfft.n_rows, input.n_cols, input.n_cols);
+  for (us i = 0; i < input.n_cols; i++) {
+    for (us j = 0; j < input.n_cols; j++) {
+      output.slice(j).col(i) = rfft.col(i) % arma::conj(rfft.col(j));
+    }
+  }
+  return output;
+}
+
+AvPowerSpectra::AvPowerSpectra(const us nfft, const Window::WindowType w,
+    const d overlap_percentage,
+    const int time_constant)
+  : _ps(nfft, w) {
+
+    DEBUGTRACE_ENTER;
+    if (overlap_percentage >= 100 || overlap_percentage < 0) {
+      throw runtime_error("Overlap percentage should be >= 0 and < 100");
+    }
+
+    _overlap_jump = nfft - (nfft * overlap_percentage) / 100;
+    if (_overlap_jump == 0) {
+      throw runtime_error("Overlap is too high. Results in no jump. Please "
+          "choose a smaller overlap percentage or a higher nfft");
+    }
+    if(time_constant < 0) {
+      _mode = Mode::Averaging;
+    } else if(time_constant == 0) {
+      _mode = Mode::Spectrogram;
+    }
+    else {
+      _mode = Mode::Leaking;
+      _alpha = exp(-static_cast<d>(_overlap_jump)/time_constant);
+    }
+
+  }
+
+arma::cx_cube* AvPowerSpectra::compute(const dmat& timedata) {
+
+
+  return nullptr;
+}

src/lasp/dsp/lasp_avpowerspectra.h

@@ -0,0 +1,155 @@
+#pragma once
+#include <memory>
+#include "lasp_fft.h"
+#include "lasp_mathtypes.h"
+#include "lasp_window.h"
+#include <optional>
+
+/**
+ * @brief Computes single-sided cross-power spectra for a group of channels
+ */
+class PowerSpectra {
+public:
+  /**
+   * @brief The FFT length
+   */
+  us nfft;
+
+private:
+  /**
+   * @brief Instance to compute actual FFT's
+   */
+  Fft _fft;
+
+  vd _window;
+  c _scale_fac;
+
+public:
+  /**
+   * @brief Initalize power spectra computer
+   *
+   * @param nfft The fft length
+   * @param w The window type
+   */
+  PowerSpectra(const us nfft, const Window::WindowType w);
+
+  /**
+   * @brief Initalize power spectra computer
+   *
+   * @param window Uses the window length to compute fft length, and uses the
+   * window shape as the actual window.
+   */
+  PowerSpectra(const vd &window);
+
+  /**
+   * @brief Computes the spectra. Data is normalized by the (spectral) power of
+   * the window, to compensate for the effect of the window.
+   *
+   * @param input Input data, first axis is assumed the sample, second axis the
+   * channel. Input first dimension should be equal to nfft, otherwise a
+   * runtime error is thrown.
+   *
+   * @return Cross-power-spectra. Cubic array with size indices C_fij, where
+   * f is the frequency index, i is the row and j is the column. An element
+   * can be accessed by: C(f, i, j). Storage is such that frequency components
+   * are contiguous.
+   */
+  arma::Cube<c> compute(const dmat &input);
+};
+
+/**
+ * @brief Estimate cross-power spectra using Welch' method of spectral
+ * estimation. The exact amount of overlap in Welch' method is rounded up to a
+ * certain amount of samples.
+ */
+class AvPowerSpectra {
+
+  enum class Mode {
+    Averaging = 0,  // Averaging all time date
+    Leaking = 1,    // Exponential weighting of an "instantaneous cps"
+    Spectrogram = 2 // Instantenous spectrum, no averaging
+  };
+
+  Mode _mode;
+  d _alpha; // Only valid in case of 'Leaking'
+
+  PowerSpectra _ps;
+  /**
+   * @brief Current estimate of cross-spectral density
+   */
+  arma::cx_cube _est;
+
+  /**
+   * @brief Buffer of storage of time data.
+   */
+  dmat _timeBuf;
+  /**
+   * @brief Current time index in ring buffer
+   */
+  us _curTimeIdx = 0;
+  /**
+   * @brief The amount of samples to jump in the overlap
+   */
+  us _overlap_jump;
+
+public:
+  /**
+   * @brief Initalize averaged power spectra computer. If a time constant is
+   * given > 0, it is used in a kind of exponential weighting.
+   *
+   * @param nfft The fft length
+   * @param w The window type.
+   * @param overlap_percentage A number 0 < overlap_percentage <= 100. It
+   * determines the amount of overlap used in Welch' method. A typical value is
+   * 50 %, i.e. 50.
+   * @param time_constant Value should either be < 0, indicating that the
+   * estimate is averages over all time data.
+   * For a value = 0 the instantaneous power spectrum is returned, which can be
+   * interpreted as the spectrogram. For a value > 0 a exponential forgetting is
+   * used, where the value is used as the time constant such that the decay
+   * follows approximately the trend exp(-n/time_constant), where n is the
+   * sample number in the power spectra. To choose 'fast' time weighting, set
+   * time_constant to the value of fs*0.125, where fs denotes the sampling
+   * frequency.
+   **/
+  AvPowerSpectra(const us nfft = 2048,
+                 const Window::WindowType w = Window::WindowType::Hann,
+                 const d overlap_percentage = 50.,
+                 const int time_constant = -1);
+
+  AvPowerSpectra(const AvPowerSpectra &) = delete;
+  AvPowerSpectra &operator=(const AvPowerSpectra &) = delete;
+
+  void reset() { _est.reset(); }
+
+  /**
+   * @brief Comnpute an update of the power spectra based on given time data.
+   * Note that the number of channels is determined from the first time this
+   * function is called. If a later call has an incompatible number of
+   * channels, a runtime error is thrown.
+   *
+   * @param timedata
+   *
+   * @return Optionally, a reference (NOT OWNING POINTER) to a new estimate of
+   * the power spectra. An update is only given if the amount of new time data
+   * is enough to compute a new estimate. This is the case when the amount of
+   * new time data is enough.
+   */
+  arma::cx_cube *compute(const dmat &timedata);
+
+  /**
+   * @brief Returns the latest estimate of cps (cross-power spectra.
+   *
+   * @return Pointer (reference, not owning!) to spectral estimate date, or
+   * nullptr, in case nothing could yet be computed.
+   */
+  arma::cx_cube *get_est();
+
+  /**
+   * @brief The overlap is rounded to a certain amount of time samples. This
+   * function returns that value.
+   *
+   * @return The amount of samples in overlapping.
+   */
+  us exactOverlapSamples() const { return _ps.nfft - _overlap_jump; }
+};

src/lasp/dsp/lasp_fft.cpp

@@ -0,0 +1,151 @@
+// lasp_fft.c
+//
+// Author: J.A. de Jong - ASCEE
+//
+// Description:
+// FFt implementation
+//
+//////////////////////////////////////////////////////////////////////
+#include <memory>
+#define DEBUGTRACE_ENABLED
+#include "lasp_fft.h"
+#include "debugtrace.hpp"
+#include "lasp_config.h"
+using std::runtime_error;
+
+#if LASP_FFT_BACKEND == Armadillo
+#include "fftpack.h"
+class Fft_impl {
+public:
+  us nfft;
+  Fft_impl(const us nfft) : nfft(nfft) {
+    throw runtime_error(
+        "This code does not output correct results, as it computes the full "
+        "FFT. It needs to be reworked to single-sided spectra.");
+
+    DEBUGTRACE_ENTER;
+  }
+
+  vd ifft(const vc &f) { return arma::ifft(f); }
+  vc fft(const vd &time) { return arma::fft(time); }
+};
+#elif LASP_FFT_BACKEND == FFTW
+#include <fftw3.h>
+
+/**
+ * @brief FFTW implementation
+ */
+class Fft_impl {
+public:
+  us nfft;
+  fftw_plan forward_plan = nullptr;
+  fftw_plan reverse_plan = nullptr;
+  fftw_complex *frequencyDomain = nullptr;
+  d *timeDomain = nullptr;
+
+  Fft_impl(const us nfft) : nfft(nfft) {
+
+    timeDomain = (d *)fftw_malloc(sizeof(d) * nfft);
+
+    frequencyDomain =
+        (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (nfft / 2 + 1));
+
+    forward_plan =
+        fftw_plan_dft_r2c_1d(nfft, timeDomain, frequencyDomain, FFTW_MEASURE);
+
+    reverse_plan =
+        fftw_plan_dft_c2r_1d(nfft, frequencyDomain, timeDomain, FFTW_MEASURE);
+  };
+
+  ~Fft_impl() {
+    fftw_destroy_plan(forward_plan);
+    fftw_destroy_plan(reverse_plan);
+    fftw_free(frequencyDomain);
+    fftw_free(timeDomain);
+  }
+
+  vc fft(const vd &time) {
+
+    vc res(nfft / 2 + 1);
+    memcpy(timeDomain, time.memptr(), sizeof(d) * nfft);
+
+    fftw_execute(forward_plan);
+
+    memcpy(reinterpret_cast<void *>(res.memptr()), frequencyDomain,
+           sizeof(d) * (nfft / 2 + 1));
+
+    return res;
+  }
+  vd ifft(const vc &f) {
+    vd res(nfft);
+    memcpy(frequencyDomain,
+           reinterpret_cast<void *>(const_cast<c *>(f.memptr())),
+           (nfft / 2 + 1) * sizeof(c));
+
+    fftw_execute(reverse_plan);
+
+    memcpy(res.memptr(), frequencyDomain, nfft * sizeof(d));
+
+    /* Scale by dividing by nfft. Checked with numpy implementation
+     * that this indeed needs to be done for FFTW. */
+    res *= 1 / nfft;
+
+    return res;
+  }
+};
+#else
+#error                                                                         \
+    "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW"
+#endif
+
+Fft::Fft(const us nfft) {
+  DEBUGTRACE_ENTER;
+  if (nfft == 0) {
+    throw std::runtime_error("Invalid nfft: 0");
+  }
+  if (nfft >= LASP_MAX_NFFT) {
+    throw std::runtime_error("Invalid nfft, should be smaller than: " +
+                             std::to_string(LASP_MAX_NFFT));
+  }
+  _impl = std::make_unique<Fft_impl>(nfft);
+}
+Fft::~Fft() { }
+
+us Fft::nfft() const { return _impl->nfft; }
+
+vd Fft::ifft(const vc &freqdata) {
+  DEBUGTRACE_ENTER;
+  if (freqdata.size() != (_impl->nfft / 2 + 1)) {
+    throw runtime_error("Invalid size for frequency domain data rows."
+                        " Should be equal to floor(nfft/2)+1");
+  }
+  return _impl->ifft(freqdata);
+}
+dmat Fft::ifft(const cmat &freqdata) {
+  dmat res(_impl->nfft, freqdata.n_cols);
+  for (us colno = 0; colno < freqdata.n_cols; colno++) {
+    res.col(colno) = _impl->ifft(freqdata.col(colno));
+  }
+  return res;
+}
+vc Fft::fft(const vd &timedata) { return _impl->fft(timedata); }
+
+void load_fft_wisdom(const char *wisdom) {
+#if LASP_FFT_BACKEND == Armadillo
+#elif LASP_FFT_BACKEND == FFTW
+  if (wisdom) {
+    int rv = fftw_import_wisdom_from_string(wisdom);
+    if (rv != 1) {
+      fprintf(stderr, "Error loading FFTW wisdom");
+    }
+  }
+#endif
+}
+
+char *store_fft_wisdom() {
+#if LASP_FFT_BACKEND == Armadillo
+  return NULL;
+#elif LASP_FFT_BACKEND == FFTW
+  return fftw_export_wisdom_to_string();
+#endif
+}

src/lasp/dsp/lasp_fft.h

@@ -0,0 +1,90 @@
+// lasp_fft.h
+//
+// Author: J.A. de Jong - ASCEE
+//
+// Description: FFT class
+// Interface to the FFT library, multiple channel FFT's
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#include <memory>
+#include "lasp_mathtypes.h"
+
+/**
+ * Load wisdom data from a NULL-terminated character string. Note that at this
+ * point in time this is only used by the FFTW fft backend.
+ *
+ * @param[in] wisdom: Wisdom string
+ */
+void load_fft_wisdom(const char* wisdom);
+
+/**
+ * Creates a NULL-terminated string containing FFTW wisdom, or state knowledge
+ * from other libraries.
+ *
+ * @returns A NULL-terminated string. *Should be deallocated after being used*
+ */
+char* store_fft_wisdom(void);
+
+class Fft_impl;
+
+/**
+ * Perform forward FFT's on real time data. Computes single-sided spectra.
+ */
+class Fft {
+  std::unique_ptr<Fft_impl> _impl;
+
+  public:
+  /**
+   * @brief Initialize FFT
+   *
+   * @param nfft The length of nfft
+   */
+    Fft(const us nfft);
+    ~Fft();
+
+    Fft(const Fft&) = delete;
+    Fft& operator=(const Fft&) = delete;
+
+    /**
+     * @brief Return nfft
+     *
+     * @return nfft
+     */
+    us nfft() const;
+
+    /** 
+     * Compute the fft for a single channel of data.
+     *
+     * @param[in] timedata Input time data, should have size nfft
+     * @returns result Result complex vector
+     * */
+    vc fft(const vd& timedata);
+
+    /** 
+     * Compute the fft of the data matrix, first axis is assumed to be
+     * the time axis.
+     *
+     * @param[in] timedata Input time data, should have size nfft. First axis
+     * is time, second axis is channel
+     * @returns Result complex array
+     */
+    cmat fft(const dmat& timedata);
+
+    /** 
+     * Perform inverse fft on a single channel.
+     * 
+     * @param[in] freqdata Frequency domain input data, to be iFft'th. Should
+     * have size nfft/2+1
+     * @returns timedata: iFft't data, size nfft.
+     */
+    vd ifft(const vc& freqdata);
+
+    /** 
+     * Perform inverse FFT
+     *
+     * @param[in] freqdata Frequency domain data
+     * @param[out] timedata Time domain result
+     */
+    dmat ifft(const cmat& freqdata);
+
+};

src/lasp/dsp/lasp_mathtypes.h

@@ -40,6 +40,8 @@
 #endif // LASP_DOUBLE_PRECISION
 
 using vd = arma::Col<d>;
+using vc = arma::Col<c>;
 using dmat = arma::Mat<d>;
+using cmat = arma::Mat<c>;
 
 const d number_pi = arma::datum::pi;

src/lasp/lasp_cpp.cpp

@@ -1,5 +1,23 @@
 #include "lasp_config.h"
 #include <pybind11/pybind11.h>
+
+/*! \mainpage
+ *
+ * \section intro_sec Introduction
+ *
+ * Welcome to the LASP (Library for Acoustic Signal Processing) code
+ * documentation. The code comprises a part which is written in C++, a part
+ * that is written in Python, and a part that functions as glue, which is
+ * Pybind11 C++ glue code. An example of such a file is the current one.
+ *
+ * \section Installation
+ *
+ * For the installation manual, please refer to the README.md of the Git
+ * repository. It is recommended to install the software from source.
+ *
+ *
+ * */
+
 namespace py = pybind11;
 
 void init_dsp(py::module &m);
@@ -18,9 +36,12 @@ PYBIND11_MODULE(lasp_cpp, m) {
   init_streammgr(m);
   init_datahandler(m);
 
+  // We store the version number of the code via CMake, and create an
+  // attribute in the C++ code.
   m.attr("__version__") = std::to_string(LASP_VERSION_MAJOR) + "." +
                           std::to_string(LASP_VERSION_MINOR);
 
+
   m.attr("LASP_VERSION_MAJOR") = LASP_VERSION_MAJOR;
   m.attr("LASP_VERSION_MINOR") = LASP_VERSION_MINOR;
 }

third_party/CMakeLists.txt

@@ -1,10 +0,0 @@
-# third_party/CMakeLists.txt
-
-if(LASP_FFT_BACKEND STREQUAL "FFTW")
-    find_library(FFTW_LIBRARY NAMES fftw3 fftw)
-    set(FFT_LIBRARIES "${FFTW_LIBRARY}")
-elseif(LASP_FFT_BACKEND STREQUAL "FFTPack")
-    include_directories(fftpack)
-    set(FFT_LIBRARIES fftpack)
-    add_subdirectory(fftpack)
-endif()

third_party/carma

@@ -0,0 +1 @@
+Subproject commit 5673bc713e9515a00fba0c7378375fc075bdc693

third_party/fftpack/CMakeLists.txt

@@ -1,9 +0,0 @@
-# We borrow Numpy's implementation for doing the Fast Fourier Transform.
-# This FFT code appears to be faster than KISSFFT.
-add_library(fftpack
-  fftpack.c
-)
-
-# Ling fft to math
-target_link_libraries(fftpack PRIVATE m)
-target_include_directories(fftpack PUBLIC .)

third_party/fftpack/fftpack.h

@@ -1,27 +0,0 @@
-/*
- * This file is part of tela the Tensor Language.
- * Copyright (c) 1994-1995 Pekka Janhunen
- */
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-#define NPY_VISIBILITY_HIDDEN
-
-#ifdef LASP_DOUBLE_PRECISION
-#define Treal double
-#else
-#define Treal float
-#endif
-
-extern NPY_VISIBILITY_HIDDEN void npy_cfftf(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_cfftb(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_cffti(int N, Treal wrk[]);
-
-extern NPY_VISIBILITY_HIDDEN void npy_rfftf(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_rfftb(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_rffti(int N, Treal wrk[]);
-
-#ifdef __cplusplus
-}
-#endif