Bump v1.6.6

.gitea/workflows/workflow.yml

@@ -12,7 +12,7 @@ jobs:
         - lasp_dist:/dist
     steps:
       - name: Checkout
-        uses: actions/checkout@v4
+        uses: actions/checkout@v3
         with:
           submodules: true
 
@@ -25,7 +25,7 @@ jobs:
 
       - name: Cleanup old dist files and copy new to /dist dir
         run: |-
-          rm /dist/*
+          rm -f /dist/*
           cp -v dist/* /dist
 
   Release-Ubuntu:

.gitignore

@@ -21,3 +21,6 @@ acme_log.log
 .venv
 .py-build-cmake_cache
 cpp_src/lasp_config.h
+.cache
+.vscode
+build

CMakeLists.txt

@@ -1,17 +1,27 @@
-cmake_minimum_required (VERSION 3.16) 
-project(LASP LANGUAGES C CXX VERSION 1.1)
+cmake_minimum_required(VERSION 3.16)
+project(LASP LANGUAGES C CXX VERSION 1.6.3)
 
 set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED)
 
-option(LASP_DOUBLE_PRECISION "Compile as double precision floating point" ON)
+if(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "aarch64")
+    set(RPI TRUE)
+else()
+    set(RPI FALSE)
+endif()
 
 # Setting defaults for PortAudio and RtAudio backend, depending on Linux /
 # Windows.
+set(DEFAULT_DOUBLE_PRECISION ON)
 if(WIN32)
     set(DEFAULT_RTAUDIO OFF)
     set(DEFAULT_PORTAUDIO ON)
     set(DEFAULT_ULDAQ OFF)
+elseif(${RPI})
+    set(DEFAULT_RTAUDIO OFF)
+    set(DEFAULT_PORTAUDIO ON)
+    set(DEFAULT_ULDAQ OFF)
+    set(DEFAULT_DOUBLE_PRECISION OFF)
 else()
     set(DEFAULT_RTAUDIO OFF)
     set(DEFAULT_PORTAUDIO ON)
@@ -19,6 +29,7 @@ else()
 endif()
 
 
+option(LASP_DOUBLE_PRECISION "Compile as double precision floating point" ${DEFAULT_DOUBLE_PRECISION})
 option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ${DEFAULT_RTAUDIO})
 option(LASP_HAS_PORTAUDIO "Compile with PortAudio Daq backend" ${DEFAULT_PORTAUDIO})
 if(LASP_HAS_PORTAUDIO AND LASP_HAS_RTAUDIO)
@@ -34,7 +45,7 @@ option(LASP_BUILD_CPP_TESTS "Build CPP test code" OFF)
 find_program(CCACHE_PROGRAM ccache)
 if(CCACHE_PROGRAM)
     set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE "${CCACHE_PROGRAM}")
-    set_property(GLOBAL PROPERTY RULE_LAUNCH_LINK  "${CCACHE_PROGRAM}")
+    set_property(GLOBAL PROPERTY RULE_LAUNCH_LINK "${CCACHE_PROGRAM}")
 endif()
 
 # To allow linking to static libs from other directories
@@ -80,7 +91,7 @@ endif()
 # Tune for current machine
 if(LASP_BUILD_TUNED)
     if(CMAKE_BUILD_TYPE STREQUAL "Debug")
-        message(FATAL_ERROR 
+        message(FATAL_ERROR
             "Cannot build optimized and tuned code when debug is switched on")
     endif()
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -mtune=native")
@@ -105,20 +116,27 @@ else()
 endif()
 
 # ###################################### Compilation flags
-set(CMAKE_C_FLAGS_RELEASE  "-O3 -flto -mfpmath=sse -march=x86-64 -mtune=native \
+set(CMAKE_C_FLAGS_RELEASE "-O3 -flto -mtune=native \
 -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions")
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-type-limits -Werror=return-type")
 
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -mtune=native \
+-fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions")
+
+if(NOT ${RPI})
+    set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} -mfpmath=sse -march=x86-64")
+    set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -mfpmath=sse -march=x86-64")
+endif()
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -Wall ")
 
 # ############################# End compilation flags
-include_directories(/usr/lib/python3.10/site-packages/numpy/core/include)
 
 # ####################################### End of user-adjustable variables section
 include(OSSpecific)
 include(rtaudio)
 include(portaudio)
 include(uldaq)
-# 
+#
 add_definitions(-Dgsl_CONFIG_DEFAULTS_VERSION=1)
 add_subdirectory(cpp_src)
 if(LASP_BUILD_CPP_TESTS)

README.md

@@ -1,4 +1,5 @@
-# Library for Acoustic Signal Processing
+Library for Acoustic Signal Processing
+======================================
 
 
 Welcome to LASP: Library for Acoustic Signal Processing. LASP is a C++ library
@@ -45,17 +46,14 @@ If you have any question(s), please feel free to contact us: [email](info@ascee.
 
 # Installation - Linux (Ubuntu-based)
 
-## From wheel (recommended for non-developers)
-
-### Prerequisites
+## Prerequisites
 
 Run the following on the command line to install all prerequisites on
 Debian-based Linux:
 
-- `sudo apt install python3-pip libfftw3-3 libopenblas-base libusb-1.0-0
-  libpulse0`
+- `sudo apt install python3-pip libfftw3-3 libopenblas-base libusb-1.0-0 libpulse0`
 
-### Download and install LASP
+## Installation from wheel (recommended for non-developers)
 
 Go to: [LASP releases](https://code.ascee.nl/ASCEE/lasp/releases/latest/) and
 download the latest `.whl`. Then run:
@@ -84,35 +82,43 @@ If building RtAudio with the Jack Audio Connection Kit (JACK) backend, you will
 - `$ cd lasp`
 - `pip install -e .`
 
-# Installation - (x86_64) Windows (with WinPython), build with MSYS2 (NOT YET UPDATED!!)
+# Installation - (x86_64) Windows (with WinPython), build with MSYS2
 
 ## Prerequisites
 
 - Download and install [WinPython](https://winpython.github.io)
+
+## From wheel
+
+- Download latest wheel from [LASP releases](https://code.ascee.nl/ASCEE/lasp/releases/latest/) and
+download the latest `.whl`. Then install with `pip`.
+
+## From source
+
 - Download and install [MSYS2](https://msys2.org). Make sure to install the
   x86_64 version.
-- Download and install [Git for Windows](https://git-scm.com)
-
 - When unzipping WinPython, make sure to choose a proper and simple path, i.e.
   C:\winpython
-- Append C:\winpython\ to the PATH environment variable.
-- Run Python and install Pybind11
-    - `python -m pip install pybind11`
-
-
-- Open a msys2 **MINGW64** terminal. And run:
-
-    - `pacman -S git`
-
-- Then clone the LASP repo:
-
-    - `git clone https://code.ascee.nl/ascee/lasp`
-    - `cd lasp`
-
-- Configure MSYS2 further, and run cmake:
-
-    - `scripts/install_msys2_buiddeps.sh`
-    - `scripts/configur_cmake_msys2.sh`
+- Download and install [Git for Windows](https://git-scm.com)
+- Open an MSYS2 **MINGW64** terminal, and install some tools we require:
+  - `$ pacman -S git`
+- Create a new virtualenv:
+  - `$ /c/winpython/<py-distr-dir>/python.exe -m venv venv`
+- Add the venv-python to the path (eases a lot of commands)
+  - `$ export PATH=$PATH:~/venv/Scripts`
+- Install `build`:
+  - `$ pip install build`
+- Clone LASP:
+  - `$ git clone --recurse-submodules https://code.ascee.nl/ascee/lasp && cd lasp`
+- If run for the first time, we have to install the libraries we depend on in
+  MSYS2 (this only has to be done on a fresh MSYS2 installation):
+  - `$ scripts/install_msys2_builddeps.sh`
+- Copy over required DLL's to be included in distribution:
+  - `scripts/copy_windows_dlls.sh`
+- And... build!
+  - `pyproject-build`
+- Lastly: the generated wheel can be installed in the current virtualenv:
+  - `pip install dist/lasp*.whl`
 
 
 # Documentation
@@ -121,7 +127,7 @@ If building RtAudio with the Jack Audio Connection Kit (JACK) backend, you will
 
 [Online LASP documentation](https://lasp.ascee.nl/).
 
-## In directory
+## In directory (Linux/Debian)
 
 `$ sudo apt install doxygen graphviz`
 `$ pip install doxypypy`
@@ -138,3 +144,25 @@ This will build the documentation. It can be read by:
 
 - See examples directories for IPython notebooks.
 - Please refer to the [documentation](https://lasp.ascee.nl/) for features.
+
+
+# Development docs
+
+## Bumping version number
+
+When bumping the version number, please update the number in
+
+- `pyproject.toml`
+- `CMakeLists.txt`
+
+Then, create a commit with tag `vX.X.X`, and push it.
+
+## Updating to latest version (editable mode)
+
+When updating to the latest version of LASP in editable mode:
+
+```bash
+- $ git pull
+- $ git submodule update
+- $ pip install -e . -v
+```

cmake/portaudio.cmake

@@ -2,6 +2,11 @@
 if(LASP_HAS_PORTAUDIO)
     message("Building with Portaudio backend")
     if(WIN32)
+        set(PA_USE_ALSA FALSE CACHE BOOL "Build PortAudio with ALSA backend")
+		set(PA_USE_ASIO TRUE CACHE BOOL "Build PortAudio with ASIO backend")
+		set(PA_USE_DS FALSE CACHE BOOL "Build PortAudio with Directsound backend")
+		set(PA_USE_WMME FALSE CACHE BOOL "Build PortAudio with WMME backend")
+		set(PA_USE_WDMKS FALSE CACHE BOOL "Build PortAudio with WDMKS backend")
     else()
         # Unix
         set(PA_USE_ALSA TRUE CACHE BOOL "Build PortAudio with ALSA backend")
@@ -9,7 +14,15 @@ if(LASP_HAS_PORTAUDIO)
         set(PA_USE_PULSEAUDIO FALSE CACHE BOOL "Build PortAudio with PulseAudio backend")
         set(PA_BUILD_SHARED_LIBS FALSE CACHE BOOL "Build static library")
     endif()
+
     add_subdirectory(third_party/portaudio)
     include_directories(third_party/portaudio/include)
     link_directories(third_party/portaudio)
+
+    if(PA_USE_ALSA)
+    add_definitions(-DLASP_HAS_PA_ALSA=1)
+    else()
+    add_definitions(-DLASP_HAS_PA_ALSA=0)
+    endif()
+
 endif()

cpp_src/CMakeLists.txt

@@ -38,18 +38,26 @@ pybind11_add_module(lasp_cpp MODULE lasp_cpp.cpp
 target_link_libraries(lasp_cpp PRIVATE lasp_device_lib lasp_dsp_lib
     ${OpenMP_CXX_LIBRARIES} ${LASP_FFT_LIBS} ${TARGET_OS_LINKLIBS})
 
+target_compile_definitions(lasp_cpp PRIVATE
+    MODULE_NAME=$<TARGET_FILE_BASE_NAME:lasp_cpp>
+    VERSION_INFO="${PY_FULL_VERSION}"
+)
+
+# Hide all symbols by default (including external libraries on Linux)
+if(CMAKE_BUILD_TYPE STREQUAL "Release")
+	set_target_properties(lasp_cpp PROPERTIES
+		CXX_VISIBILITY_PRESET "hidden"
+		VISIBILITY_INLINES_HIDDEN true)
+	if (CMAKE_SYSTEM_NAME MATCHES "Linux")
+		target_link_options(lasp_cpp PRIVATE "LINKER:--exclude-libs,ALL")
+	endif()
+endif()
+
 if(DEFINED PY_BUILD_CMAKE_MODULE_NAME)
     # Install the Python module
     install(TARGETS lasp_cpp
             EXCLUDE_FROM_ALL
             COMPONENT python_modules
             DESTINATION ${PY_BUILD_CMAKE_MODULE_NAME})
-    # Install the debug file for the Python module (Windows only)
-    if (WIN32)
-        install(FILES $<TARGET_PDB_FILE:_add_module>
-                EXCLUDE_FROM_ALL
-                COMPONENT python_modules
-                DESTINATION ${PY_BUILD_CMAKE_MODULE_NAME}
-                OPTIONAL)
-    endif()
+
 endif()

cpp_src/device/lasp_daq.cpp

@@ -1,4 +1,4 @@
-/* #define DEBUGTRACE_ENABLED */
+// #define DEBUGTRACE_ENABLED
 #include "debugtrace.hpp"
 #include "lasp_daqconfig.h"
 

cpp_src/device/lasp_indatahandler.cpp

@@ -1,11 +1,11 @@
-/* #define DEBUGTRACE_ENABLED */
+// #define DEBUGTRACE_ENABLED
 #include "lasp_indatahandler.h"
 #include "debugtrace.hpp"
 #include "lasp_streammgr.h"
 #include <thread>
 
 InDataHandler::InDataHandler(SmgrHandle mgr, const InCallbackType cb,
-                             const InResetType resetfcn)
+                             const ResetCallbackType resetfcn)
     : _mgr(mgr), inCallback(cb), reset(resetfcn)
 #if LASP_DEBUG == 1
       ,
@@ -29,20 +29,22 @@ void InDataHandler::start() {
 }
 void InDataHandler::stop() {
   DEBUGTRACE_ENTER;
-  checkRightThread();
+  // checkRightThread();
 #if LASP_DEBUG == 1
   stopCalled = true;
 #endif
-  if (SmgrHandle handle = _mgr.lock()) {
-    handle->removeInDataHandler(*this);
+  if (SmgrHandle smgr = _mgr.lock()) {
+    smgr->removeInDataHandler(*this);
+  } else {
+    DEBUGTRACE_PRINT("No stream manager alive anymore!");
   }
 }
 
 InDataHandler::~InDataHandler() {
 
   DEBUGTRACE_ENTER;
-  checkRightThread();
 #if LASP_DEBUG == 1
+  // checkRightThread();
   if (!stopCalled) {
     std::cerr << "************ BUG: Stop function not called while arriving at "
                  "InDataHandler's destructor. Fix this by calling "

cpp_src/device/lasp_indatahandler.h

@@ -22,7 +22,7 @@ using InCallbackType = std::function<void(const DaqData &)>;
 /**
  * @brief Function definition for the reset callback.
  */
-using InResetType = std::function<void(const Daq *)>;
+using ResetCallbackType = std::function<void(const Daq *)>;
 
 class InDataHandler {
 
@@ -38,7 +38,7 @@ protected:
 public:
   ~InDataHandler();
   const InCallbackType inCallback;
-  const InResetType reset;
+  const ResetCallbackType reset;
 
   /**
    * @brief When constructed, the handler is added to the stream manager, which
@@ -50,7 +50,7 @@ public:
    * changes state.
    */
   InDataHandler(SmgrHandle mgr, InCallbackType cb,
-                InResetType resetfcn);
+                ResetCallbackType resetfcn);
 
   /**
    * @brief Adds the current InDataHandler to the list of handlers in the

cpp_src/device/lasp_streammgr.cpp

@@ -1,4 +1,4 @@
-/* #define DEBUGTRACE_ENABLED */
+// #define DEBUGTRACE_ENABLED
 #include "lasp_streammgr.h"
 
 #include <assert.h>
@@ -8,12 +8,15 @@
 #include <iostream>
 #include <memory>
 #include <mutex>
+#include <thread>
 
 #include "debugtrace.hpp"
 #include "lasp_biquadbank.h"
 #include "lasp_indatahandler.h"
 #include "lasp_thread.h"
 
+using namespace std::literals::chrono_literals;
+
 using std::cerr;
 using std::endl;
 using rte = std::runtime_error;
@@ -27,7 +30,7 @@ using rte = std::runtime_error;
 std::weak_ptr<StreamMgr> _mgr;
 std::mutex _mgr_mutex;
 
-using Lck = std::scoped_lock<std::mutex>;
+using Lck = std::scoped_lock<std::recursive_mutex>;
 
 /**
  * @brief The only way to obtain a stream manager, can only be called from the
@@ -38,11 +41,11 @@ using Lck = std::scoped_lock<std::mutex>;
 SmgrHandle StreamMgr::getInstance() {
   DEBUGTRACE_ENTER;
 
+  std::scoped_lock<std::mutex> lck(_mgr_mutex);
   auto mgr = _mgr.lock();
   if (!mgr) {
     // Double Check Locking Pattern, if two threads would simultaneously
     // instantiate the singleton instance.
-    Lck lck(_mgr_mutex);
 
     auto mgr = _mgr.lock();
     if (mgr) {
@@ -84,37 +87,70 @@ void StreamMgr::rescanDAQDevices(bool background,
                                  std::function<void()> callback) {
   DEBUGTRACE_ENTER;
   DEBUGTRACE_PRINT(background);
+  if (_scanningDevices) {
+    throw rte("A background device scan is already busy");
+  }
 
+  Lck lck(_mtx);
   checkRightThread();
   if (_inputStream || _outputStream) {
     throw rte("Rescanning DAQ devices only possible when no stream is running");
   }
-  if (!_devices_mtx.try_lock()) {
-    throw rte("A background DAQ device scan is probably already running");
-  }
-  _devices_mtx.unlock();
 
-  std::scoped_lock lck(_devices_mtx);
   _devices.clear();
   if (!background) {
+    _scanningDevices = true;
     rescanDAQDevices_impl(callback);
   } else {
     DEBUGTRACE_PRINT("Rescanning DAQ devices on different thread...");
+    _scanningDevices = true;
     _pool.push_task(&StreamMgr::rescanDAQDevices_impl, this, callback);
   }
 }
+#if LASP_HAS_PORTAUDIO && LASP_HAS_PA_ALSA
+#include <alsa/asoundlib.h>
+void empty_handler(const char *file, int line, const char *function, int err,
+                   const char *fmt, ...) {}
+
+// Temporarily set the ALSA eror handler to something that does nothing, to
+// prevent ALSA from spitting out all kinds of misconfiguration errors.
+class MuteErrHandler {
+ private:
+  snd_lib_error_handler_t _default_handler;
+
+ public:
+  explicit MuteErrHandler() {
+    _default_handler = snd_lib_error;
+    snd_lib_error_set_handler(empty_handler);
+  }
+
+  ~MuteErrHandler() { snd_lib_error_set_handler(_default_handler); }
+};
+#else
+// Does nothin in case of no ALSA
+class MuteErrHandler {};
+#endif
+
 void StreamMgr::rescanDAQDevices_impl(std::function<void()> callback) {
   DEBUGTRACE_ENTER;
-  std::scoped_lock lck(_devices_mtx);
-  _devices = DeviceInfo::getDeviceInfo();
+  assert(!_inputStream && !_outputStream);
+  Lck lck(_mtx);
+  // Alsa spits out annoying messages that are not useful
+  {
+    MuteErrHandler guard;
+
+    _devices = DeviceInfo::getDeviceInfo();
+  }
   if (callback) {
     callback();
   }
+
+  _scanningDevices = false;
 }
 void StreamMgr::inCallback(const DaqData &data) {
   DEBUGTRACE_ENTER;
 
-  std::scoped_lock<std::mutex> lck(_inDataHandler_mtx);
+  Lck lck(_mtx);
 
   assert(_inputFilters.size() == data.nchannels);
 
@@ -139,12 +175,13 @@ void StreamMgr::inCallback(const DaqData &data) {
       }
     }
 
+    DEBUGTRACE_PRINT("Calling incallback for handlers (filtered)...");
     for (auto &handler : _inDataHandlers) {
       handler->inCallback(input_filtered);
     }
-
   } else {
     /// No input filters
+    DEBUGTRACE_PRINT("Calling incallback for handlers...");
     for (auto &handler : _inDataHandlers) {
       handler->inCallback(data);
     }
@@ -155,8 +192,7 @@ void StreamMgr::setSiggen(std::shared_ptr<Siggen> siggen) {
   DEBUGTRACE_ENTER;
   checkRightThread();
 
-  std::scoped_lock<std::mutex> lck(_siggen_mtx);
-
+  Lck lck(_mtx);
   // If not set to nullptr, and a stream is running, we update the signal
   // generator by resetting it.
   if (isStreamRunningOK(StreamType::output) && siggen) {
@@ -213,9 +249,9 @@ bool fillData(DaqData &data, const vd &signal) {
   return true;
 }
 void StreamMgr::outCallback(DaqData &data) {
-  /* DEBUGTRACE_ENTER; */
+  DEBUGTRACE_ENTER;
 
-  std::scoped_lock<std::mutex> lck(_siggen_mtx);
+  Lck lck(_mtx);
 
   if (_siggen) {
     vd signal = _siggen->genSignal(data.nframes);
@@ -244,7 +280,17 @@ void StreamMgr::outCallback(DaqData &data) {
 
 StreamMgr::~StreamMgr() {
   DEBUGTRACE_ENTER;
-  checkRightThread();
+  while (_scanningDevices) {
+    std::this_thread::sleep_for(10us);
+  }
+
+#if LASP_DEBUG == 1
+  {  // Careful, this lock needs to be released to make sure the streams can
+    // obtain a lock to the stream manager.
+    Lck lck(_mtx);
+    checkRightThread();
+  }
+#endif
   // Stream manager now handled by shared pointer. Each indata handler gets a
   // shared pointer to the stream manager, and stores a weak pointer to it.
   // Hence, we do not have to do any cleanup here. It also makes sure that the
@@ -260,13 +306,21 @@ StreamMgr::~StreamMgr() {
 }
 void StreamMgr::stopAllStreams() {
   DEBUGTRACE_ENTER;
-  checkRightThread();
+  {
+    Lck lck(_mtx);
+    checkRightThread();
+  }
+  // No lock here!
   _inputStream.reset();
   _outputStream.reset();
 }
 
 void StreamMgr::startStream(const DaqConfiguration &config) {
   DEBUGTRACE_ENTER;
+  if (_scanningDevices) {
+    throw rte("DAQ device scan is busy. Cannot start stream.");
+  }
+  Lck lck(_mtx);
   checkRightThread();
 
   bool isInput = std::count_if(config.inchannel_config.cbegin(),
@@ -278,8 +332,6 @@ void StreamMgr::startStream(const DaqConfiguration &config) {
                                 [](auto &i) { return i.enabled; }) > 0;
 
   // Find the first device that matches with the configuration
-  std::scoped_lock lck(_devices_mtx);
-
   DeviceInfo *devinfo = nullptr;
 
   // Match configuration to a device in the list of devices
@@ -404,37 +456,49 @@ void StreamMgr::startStream(const DaqConfiguration &config) {
 void StreamMgr::stopStream(const StreamType t) {
   DEBUGTRACE_ENTER;
   checkRightThread();
+  bool resetHandlers = false;
+  std::unique_ptr<Daq> *streamToStop = nullptr;
 
-  if (t == StreamType::input) {
-    if (!_inputStream) {
-      throw rte("Input stream is not running");
+  {  // Mutex locked in this scope
+    Lck lck(_mtx);
+    if (t == StreamType::input) {
+      if (!_inputStream) {
+        throw rte("Input stream is not running");
+      }
+      streamToStop = std::addressof(_inputStream);
+      resetHandlers = true;
+    } else {
+      /// t == output
+      /// Kill input stream in case that one is a duplex stream
+      if (_inputStream && _inputStream->duplexMode()) {
+        streamToStop = std::addressof(_inputStream);
+      } else {
+        if (!_outputStream) {
+          throw rte("Output stream is not running");
+        }
+        streamToStop = std::addressof(_outputStream);
+      }  // end else
     }
-    /// Kills input stream
-    _inputStream.reset();
-    /// Send reset to all in data handlers
+  }  // End of mutex lock. When stopping stream, mutex should be unlocked.
+
+  // If we arrive here, we should have a stream to stop.
+  assert(streamToStop != nullptr);
+  streamToStop->reset();
+
+  /// Send reset to all in data handlers
+  if (resetHandlers) {
+    Lck lck(_mtx);
     for (auto &handler : _inDataHandlers) {
       handler->reset(nullptr);
     }
-  } else {
-    /// t == output
-
-    /// Kill input stream in case that one is a duplex stream
-    if (_inputStream && _inputStream->duplexMode()) {
-      _inputStream.reset();
-    } else {
-      if (!_outputStream) {
-        throw rte("Output stream is not running");
-      }
-      _outputStream.reset();
-    }  // end else
   }
 }
 
 void StreamMgr::addInDataHandler(InDataHandler *handler) {
   DEBUGTRACE_ENTER;
+  Lck lck(_mtx);
   checkRightThread();
   assert(handler);
-  std::scoped_lock<std::mutex> lck(_inDataHandler_mtx);
   handler->reset(_inputStream.get());
 
   if (std::find(_inDataHandlers.cbegin(), _inDataHandlers.cend(), handler) !=
@@ -449,16 +513,17 @@ void StreamMgr::addInDataHandler(InDataHandler *handler) {
 
 void StreamMgr::removeInDataHandler(InDataHandler &handler) {
   DEBUGTRACE_ENTER;
-  checkRightThread();
-  std::scoped_lock<std::mutex> lck(_inDataHandler_mtx);
+  Lck lck(_mtx);
+  // checkRightThread();
   _inDataHandlers.remove(&handler);
 
   DEBUGTRACE_PRINT(_inDataHandlers.size());
 }
 
 Daq::StreamStatus StreamMgr::getStreamStatus(const StreamType type) const {
-  /* DEBUGTRACE_ENTER; */
+  DEBUGTRACE_ENTER;
 
+  Lck lck(_mtx);
   checkRightThread();
   // Default constructor, says stream is not running, but also no errors
 
@@ -471,6 +536,7 @@ Daq::StreamStatus StreamMgr::getStreamStatus(const StreamType type) const {
 }
 
 const Daq *StreamMgr::getDaq(StreamType type) const {
+  Lck lck(_mtx);
   checkRightThread();
 
   if (type == StreamType::input) {

cpp_src/device/lasp_streammgr.h

@@ -1,19 +1,19 @@
 #pragma once
-#include "lasp_daq.h"
-#include "lasp_siggen.h"
-#include "lasp_thread.h"
 #include <list>
 #include <memory>
 #include <mutex>
 #include <thread>
 
+#include "lasp_daq.h"
+#include "lasp_siggen.h"
+#include "lasp_thread.h"
+
 /** \addtogroup device
  * @{
  */
 class StreamMgr;
 class InDataHandler;
 
-
 class SeriesBiquad;
 
 /**
@@ -25,12 +25,15 @@ class SeriesBiquad;
  * fact is asserted.
  */
 class StreamMgr {
+  mutable std::recursive_mutex _mtx;
 
   /**
    * @brief Storage for streams.
    */
   std::unique_ptr<Daq> _inputStream, _outputStream;
 
+  std::atomic<bool> _scanningDevices{false};
+
   GlobalThreadPool _pool;
 
   /**
@@ -39,22 +42,18 @@ class StreamMgr {
    * thread-safety.
    */
   std::list<InDataHandler *> _inDataHandlers;
-  mutable std::mutex _inDataHandler_mtx;
 
   /**
    * @brief Signal generator in use to generate output data. Currently
    * implemented as to generate the same data for all output channels.
    */
   std::shared_ptr<Siggen> _siggen;
-  std::mutex _siggen_mtx;
 
   /**
    * @brief Filters on input stream. For example, a digital high pass filter.
    */
   std::vector<std::unique_ptr<SeriesBiquad>> _inputFilters;
 
-
-  mutable std::recursive_mutex _devices_mtx;
   /**
    * @brief Current storage for the device list
    */
@@ -67,9 +66,7 @@ class StreamMgr {
   friend class InDataHandler;
   friend class Siggen;
 
-
-  public:
-
+ public:
   ~StreamMgr();
 
   enum class StreamType : us {
@@ -100,9 +97,10 @@ class StreamMgr {
    * @return A copy of the internal stored list of devices
    */
   DeviceInfoList getDeviceInfo() const {
-  std::scoped_lock lck(_devices_mtx);
+    std::scoped_lock lck(_mtx);
     DeviceInfoList d2;
-    for(const auto& dev: _devices) {
+    for (const auto &dev : _devices) {
+      assert(dev != nullptr);
       d2.push_back(dev->clone());
     }
     return d2;
@@ -118,9 +116,9 @@ class StreamMgr {
    * set to true, the function returns immediately.
    * @param callback Function to call when complete.
    */
-  void
-  rescanDAQDevices(bool background = false,
-                   std::function<void()> callback = std::function<void()>());
+  void rescanDAQDevices(
+      bool background = false,
+      std::function<void()> callback = std::function<void()>());
 
   /**
    * @brief Start a stream based on given configuration.
@@ -141,12 +139,12 @@ class StreamMgr {
   }
   bool isStreamRunning(const StreamType type) const {
     switch (type) {
-    case (StreamType::input):
-      return bool(_inputStream);
-      break;
-    case (StreamType::output):
-      return bool(_outputStream);
-      break;
+      case (StreamType::input):
+        return bool(_inputStream);
+        break;
+      case (StreamType::output):
+        return bool(_outputStream);
+        break;
     }
     return false;
   }
@@ -193,11 +191,10 @@ class StreamMgr {
    */
   void setSiggen(std::shared_ptr<Siggen> s);
 
-private:
+ private:
   void inCallback(const DaqData &data);
   void outCallback(DaqData &data);
 
-
   /**
    * @brief Add an input data handler. The handler's inCallback() function is
    * called with data when available. This function should *NOT* be called by

cpp_src/device/portaudio/lasp_portaudiodaq.cpp

@@ -1,4 +1,4 @@
-/* #define DEBUGTRACE_ENABLED */
+// #define DEBUGTRACE_ENABLED
 #include "debugtrace.hpp"
 #include "lasp_config.h"
 
@@ -93,6 +93,7 @@ void fillPortAudioDeviceInfo(DeviceInfoList &devinfolist) {
           d.api = portaudioALSAApi;
           break;
         case paASIO:
+          hasDuplexMode = true;
           d.api = portaudioASIOApi;
           break;
         case paDirectSound:
@@ -422,6 +423,7 @@ Daq::StreamStatus PortAudioDaq::getStreamStatus() const {
 }
 
 PortAudioDaq::~PortAudioDaq() {
+  DEBUGTRACE_ENTER;
   PaError err;
   assert(_stream);
   if (Pa_IsStreamActive(_stream)) {
@@ -445,7 +447,7 @@ int PortAudioDaq::memberPaCallback(const void *inputBuffer, void *outputBuffer,
                                    unsigned long framesPerBuffer,
                                    const PaStreamCallbackTimeInfo *timeInfo,
                                    PaStreamCallbackFlags statusFlags) {
-  // DEBUGTRACE_ENTER;
+  DEBUGTRACE_ENTER;
   typedef Daq::StreamStatus::StreamError se;
   if (statusFlags & paPrimingOutput) {
     // Initial output buffers generated. So nothing with input yet

cpp_src/dsp/CMakeLists.txt

@@ -16,6 +16,7 @@ set(lasp_dsp_files
     lasp_threadedindatahandler.cpp
     lasp_ppm.cpp
     lasp_clip.cpp
+    lasp_freqsmooth.cpp
     )
 
 

cpp_src/dsp/lasp_avpowerspectra.cpp

@@ -2,12 +2,9 @@
 #include "lasp_avpowerspectra.h"
 #include "debugtrace.hpp"
 #include "lasp_mathtypes.h"
-#include <cmath>
-#include <optional>
+#include <stdexcept>
 
 using rte = std::runtime_error;
-using std::cerr;
-using std::endl;
 
 PowerSpectra::PowerSpectra(const us nfft, const Window::WindowType w)
     : PowerSpectra(Window::create(w, nfft)) {}

cpp_src/dsp/lasp_avpowerspectra.h

@@ -3,8 +3,6 @@
 #include "lasp_mathtypes.h"
 #include "lasp_timebuffer.h"
 #include "lasp_window.h"
-#include <memory>
-#include <optional>
 
 /** \defgroup dsp Digital Signal Processing utilities
  * These are classes and functions used for processing raw signal data, to

cpp_src/dsp/lasp_biquadbank.cpp

@@ -2,7 +2,6 @@
 #include "lasp_biquadbank.h"
 #include "debugtrace.hpp"
 #include "lasp_thread.h"
-#include <cassert>
 #include <vector>
 
 using std::cerr;

cpp_src/dsp/lasp_freqsmooth.cpp

@@ -0,0 +1,133 @@
+// #define DEBUGTRACE_ENABLED
+#include "lasp_freqsmooth.h"
+
+#include <cassert>
+
+#include "debugtrace.hpp"
+
+using rte = std::runtime_error;
+
+vd freqSmooth(const vd& freq, const vd& X, const unsigned w,
+              bool power_correct) {
+  DEBUGTRACE_ENTER;
+  if (freq.size() < 2) {
+    throw rte("Invalid frequency size. Should be > 2");
+  }
+  if (freq.size() != X.size()) {
+    throw rte("Sizes of freq and X do not match");
+  }
+  if (freq.size() > std::numeric_limits<long>::max() / 2) {
+    throw rte("Frequency size limit for smoothing is 2^30");
+  }
+  if (w == 0) {
+    throw rte("Invalid number of octaves");
+  }
+  const us Nfreq = freq.size();
+
+  // Smoothing width in unit of number of octaves
+  const d Delta = 1 / d(w);
+
+  // Minimum frequency and maximum frequency to smooth on (frequency range that
+  // is interpolated to a log scale)
+  d freq_min;
+  const d freq_max = freq(Nfreq - 1);
+  const bool firstFreqEqZero = (d_abs(freq(0)) < 1e-15);
+
+  // AC-signal power
+  d ac_pwr;
+  if (firstFreqEqZero) {
+    freq_min = freq(1);
+    if (power_correct) {
+      ac_pwr = arma::sum(X.subvec(1, Nfreq - 1));
+    }
+  } else {
+    freq_min = freq(0);
+    if (power_correct) {
+      ac_pwr = arma::sum(X);
+    }
+  }
+  DEBUGTRACE_PRINT(freq_min);
+  DEBUGTRACE_PRINT(freq_max);
+  const vd freq_log =
+      arma::logspace(d_log10(freq_min), d_log10(freq_max), 10 * Nfreq);
+  DEBUGTRACE_PRINT("freq_log = ");
+
+  const long Nfreq_sm = freq_log.size();
+
+  // Interpolate X to logscale
+  vd X_log;
+  DEBUGTRACE_PRINT("X_log = :");
+  arma::interp1(freq, X, freq_log, X_log, "*linear");
+
+  // First and last point are not interpolated well, could be minimally out of
+  // the interpolation range, due to roundoff errors. Armadillo sets these
+  // points to nan, so we have to manually "interpolate" them.
+  X_log(Nfreq_sm - 1) = X(X.size() - 1);
+  if (firstFreqEqZero) {
+    X_log(0) = X(1);
+  } else {
+    X_log(0) = X(0);
+  }
+
+  // Allocate space for smoothed X on log scale
+  vd Xsm_log(freq_log.size());
+  const d beta = d_log10(Nfreq_sm) / d_log10(2) / (Nfreq_sm - 1);
+  // int rounds down
+  const long mu = int(Delta / d(2) / beta);
+  DEBUGTRACE_PRINT(mu);
+
+  // Long is at least 32 bits. So +/- 2M points length
+  for (long k = 0; k < Nfreq_sm; k++) {
+    // const d fcur = freq_log(k);
+    long idx_start = std::max(k - mu, 0l);
+    long idx_stop = std::min(k + mu, Nfreq_sm - 1);
+
+    // Make window smaller at the sides (close to the end of the array)
+    if (idx_start == 0 || idx_stop == Nfreq_sm - 1) {
+      const long mu_edge = std::min(k - idx_start, idx_stop - k);
+      idx_start = k - mu_edge;
+      idx_stop = k + mu_edge;
+    }
+    assert(idx_stop < Nfreq_sm);
+    assert(idx_start < Nfreq_sm);
+
+    DEBUGTRACE_PRINT(idx_start)
+    DEBUGTRACE_PRINT(idx_stop);
+
+    Xsm_log(k) = arma::mean(X_log.subvec(idx_start, idx_stop));
+  }
+  DEBUGTRACE_PRINT("Xsm_log:");
+  //   std::cerr << Xsm_log << std::endl;
+
+  // Back-interpolate to a linear scale, and be wary of nans at the start end
+  // and range. Also interpolates power
+  vd Xsm(Nfreq);
+  if (firstFreqEqZero) {
+    vd Xsm_gt0;
+    arma::interp1(freq_log, Xsm_log, freq.subvec(1, Nfreq - 1), Xsm_gt0,
+                  "*linear");
+    Xsm(0) = X(0);
+    Xsm.subvec(1, Nfreq - 1) = Xsm_gt0;
+    Xsm(1) = Xsm_log(1);
+    Xsm(Nfreq - 1) = Xsm_log(Nfreq_sm - 1);
+
+    // Final step: power-correct smoothed spectrum
+    if (power_correct) {
+      d new_acpwr = arma::sum(Xsm.subvec(1, Nfreq - 1));
+      Xsm.subvec(1, Nfreq - 1) *= ac_pwr / new_acpwr;
+    }
+
+  } else {
+    arma::interp1(freq_log, Xsm_log, freq, Xsm, "*linear");
+    Xsm(0) = X(0);
+    Xsm(Nfreq - 1) = Xsm_log(Nfreq_sm - 1);
+
+    // Final step: power-correct smoothed spectrum
+    if (power_correct) {
+      d new_acpwr = arma::sum(Xsm);
+      Xsm *= ac_pwr / new_acpwr;
+    }
+  }
+
+  return Xsm;
+}

cpp_src/dsp/lasp_freqsmooth.h

@@ -0,0 +1,28 @@
+#pragma once
+#include <memory>
+#include <vector>
+
+#include "lasp_mathtypes.h"
+#include "lasp_types.h"
+
+/**
+ * \addtogroup dsp
+ * @{
+ */
+
+/**
+ * @brief Apply frequency domain smoothing to a Frequency domain (single
+ * sided)signal power spectrum
+ * 
+ * @param freq Frequency range
+ * @param X Signal pwr
+ * @param w Parameter determining the smoothing with. 1 = 1/1 octave, 3 = 1/3th
+ * octave and so on
+ * @param power_correct Apply a correction to the whole spectrum to make the
+ * signal power equal to the unsmoothed signal power.
+ * @return vd Smoothed spectrum
+ */
+vd freqSmooth(const vd& freq, const vd& X, const unsigned w,
+              bool power_correct = false);
+
+/** @} */

cpp_src/dsp/lasp_rtaps.cpp

@@ -37,7 +37,7 @@ void RtAps::inCallback(const DaqData &data) {
     cerr << "**** Error: sensitivity size does not match! *****" << endl;
     return;
   }
-  fltdata.each_row() %= _sens.as_row();
+  fltdata.each_row() /= _sens.as_row();
 
   if (_filterPrototype) {
 

cpp_src/dsp/lasp_rtsignalviewer.h

@@ -4,12 +4,9 @@
 //
 // Description: Real Time Signal Viewer.
 #pragma once
-#include "lasp_avpowerspectra.h"
-#include "lasp_filter.h"
 #include "lasp_mathtypes.h"
 #include "lasp_threadedindatahandler.h"
 #include "lasp_timebuffer.h"
-#include <memory>
 #include <mutex>
 
 /**

cpp_src/dsp/lasp_siggen.cpp

@@ -10,12 +10,12 @@ using rte = std::runtime_error;
 
 inline d level_amp(d level_dB) { return pow(10, level_dB / 20); }
 
-using mutexlock = std::scoped_lock<std::mutex>;
+using slock = std::scoped_lock<std::recursive_mutex>;
 
 vd Siggen::genSignal(const us nframes) {
 
   DEBUGTRACE_ENTER;
-  mutexlock lck(_mtx);
+  slock lck(_mtx);
 
   DEBUGTRACE_PRINT(nframes);
   vd signal(nframes, arma::fill::value(_dc_offset));
@@ -52,7 +52,7 @@ vd Siggen::genSignal(const us nframes) {
   return signal;
 }
 void Siggen::setInterruptPeriod(const d newPeriod) {
-  mutexlock lck(_mtx);
+  slock lck(_mtx);
   if (newPeriod == 0) {
     throw rte("Interruption period cannot be 0");
   }
@@ -65,7 +65,7 @@ void Siggen::setInterruptPeriod(const d newPeriod) {
 void Siggen::setFilter(const std::string &name,
                        std::shared_ptr<Filter> filter) {
   DEBUGTRACE_ENTER;
-  mutexlock lck(_mtx);
+  slock lck(_mtx);
   if (filter) {
     _filters[name] = filter;
   } else if (_filters.find(name) != _filters.end()) {
@@ -74,17 +74,17 @@ void Siggen::setFilter(const std::string &name,
 }
 void Siggen::setDCOffset(const d offset) {
   DEBUGTRACE_ENTER;
-  mutexlock lck(_mtx);
+  slock lck(_mtx);
   _dc_offset = offset;
 }
 void Siggen::setLevel(const d level, bool dB) {
   DEBUGTRACE_ENTER;
-  mutexlock lck(_mtx);
+  slock lck(_mtx);
   _level_linear = dB ? level_amp(level) : level;
 }
 void Siggen::reset(const d newFs) {
   DEBUGTRACE_ENTER;
-  mutexlock lck(_mtx);
+  slock lck(_mtx);
   _fs = newFs;
   for (auto &f : _filters) {
     assert(f.second);

cpp_src/dsp/lasp_siggen.h

@@ -27,7 +27,7 @@ private:
   bool _muted = false;
 
 protected:
-  std::mutex _mtx;
+  mutable std::recursive_mutex _mtx;
   d _fs = 0;
   /**
    * @brief Interuption of period the signal. If set, the signal will be

cpp_src/dsp/lasp_siggen_impl.cpp

@@ -7,11 +7,14 @@
 //////////////////////////////////////////////////////////////////////
 /* #define DEBUGTRACE_ENABLED */
 #include "lasp_siggen_impl.h"
-#include "debugtrace.hpp"
-#include "lasp_mathtypes.h"
+
 #include <cassert>
 
+#include "debugtrace.hpp"
+#include "lasp_mathtypes.h"
+
 using rte = std::runtime_error;
+using slock = std::scoped_lock<std::recursive_mutex>;
 
 DEBUGTRACE_VARIABLES;
 
@@ -30,6 +33,7 @@ Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) { DEBUGTRACE_ENTER; }
 
 vd Sine::genSignalUnscaled(const us nframes) {
   /* DEBUGTRACE_ENTER; */
+  slock lck(_mtx);
   const d pi = arma::datum::pi;
   vd phase_vec =
       arma::linspace(phase, phase + omg * (nframes - 1) / _fs, nframes);
@@ -41,8 +45,8 @@ vd Sine::genSignalUnscaled(const us nframes) {
 }
 
 vd Periodic::genSignalUnscaled(const us nframes) {
-
   vd res(nframes);
+  slock lck(_mtx);
   if (_signal.size() == 0) {
     throw rte("No signal defined while calling");
   }
@@ -74,15 +78,15 @@ Sweep::Sweep(const d fl, const d fu, const d Ts, const d Tq, const us flags)
 }
 
 void Sweep::resetImpl() {
-
   DEBUGTRACE_ENTER;
+  slock lck(_mtx);
 
   _cur_pos = 0;
 
   bool forward_sweep = flags & ForwardSweep;
   bool backward_sweep = flags & BackwardSweep;
 
-  const d Dt = 1 / _fs; // Deltat
+  const d Dt = 1 / _fs;  // Deltat
 
   // Estimate N, the number of samples in the sweep part (non-quiescent part):
   const us Ns = (us)(Ts * _fs);
@@ -166,7 +170,6 @@ void Sweep::resetImpl() {
       /* dVARTRACE(15, phase); */
     }
   } else if (flags & LogSweep) {
-
     DEBUGTRACE_PRINT("Log sweep");
     if (forward_sweep || backward_sweep) {
       /* Forward or backward sweep */
@@ -194,7 +197,6 @@ void Sweep::resetImpl() {
         phase += 2 * number_pi * Dt * fn;
       }
     } else {
-
       DEBUGTRACE_PRINT("Continuous sweep");
 
       const us Nf = Ns / 2;
@@ -249,17 +251,15 @@ void Sweep::resetImpl() {
         /* dbgassert(fn >= 0, "BUG"); */
 
         phase += 2 * number_pi * Dt * fn;
-        while (phase > 2 * number_pi)
-          phase -= 2 * number_pi;
+        while (phase > 2 * number_pi) phase -= 2 * number_pi;
         /* dVARTRACE(17, phase); */
       }
       /* This should be a very small number!! */
       DEBUGTRACE_PRINT(phase);
     }
-  } // End of log sweep
+  }  // End of log sweep
   else {
     // Either log or linear sweep had to be given as flags.
     assert(false);
-
   }
 }

cpp_src/dsp/lasp_siggen_impl.h

@@ -7,7 +7,7 @@
  */
 
 /**
- * \addtogroup siggen 
+ * \addtogroup siggen
  * @{
  */
 
@@ -18,8 +18,8 @@ class Noise : public Siggen {
   d level_linear;
   virtual vd genSignalUnscaled(const us nframes) override;
   void resetImpl() override;
-  public:
 
+ public:
   /**
    * @brief Constructs a noise generator. If no filter is used, the output will
    * be white noise. By default, the output will be standard deviation = 1
@@ -28,7 +28,6 @@ class Noise : public Siggen {
    */
   Noise();
   ~Noise() = default;
-
 };
 
 /**
@@ -37,13 +36,12 @@ class Noise : public Siggen {
 class Sine : public Siggen {
   d phase = 0;
   d omg;
-  protected:
 
-  void resetImpl() override final { phase=0; }
+ protected:
+  void resetImpl() override final { phase = 0; }
   virtual vd genSignalUnscaled(const us nframes) override final;
 
-  public:
-
+ public:
   /**
    * @brief Create a sine wave generator
    *
@@ -51,7 +49,7 @@ class Sine : public Siggen {
    */
   Sine(const d freq_Hz);
   ~Sine() = default;
-  void setFreq(const d newFreq) { omg = 2*arma::datum::pi*newFreq; } ;
+  void setFreq(const d newFreq) { omg = 2 * arma::datum::pi * newFreq; };
 };
 
 /**
@@ -86,8 +84,7 @@ class Sweep : public Periodic {
 
   void resetImpl() override;
 
-  public:
-
+ public:
   static constexpr int ForwardSweep = 1 << 0;
   static constexpr int BackwardSweep = 1 << 1;
   static constexpr int LinearSweep = 1 << 2;
@@ -103,11 +100,11 @@ class Sweep : public Periodic {
    * avoid temporal aliasing in case of measuring impulse responses.
    * @param[in] sweep_flags: Sweep period [s]
    */
-  Sweep(const d fl, const d fu, const d Ts, const d Tq,
-      const us sweep_flags);
+  Sweep(const d fl, const d fu, const d Ts, const d Tq, const us sweep_flags);
 
   ~Sweep() = default;
 
+
 };
 /** @} */
 /** @} */

cpp_src/dsp/lasp_threadedindatahandler.cpp

@@ -1,13 +1,15 @@
-/* #define DEBUGTRACE_ENABLED */
+// #define DEBUGTRACE_ENABLED
 #include "lasp_threadedindatahandler.h"
-#include "debugtrace.hpp"
-#include "lasp_daqdata.h"
-#include "lasp_thread.h"
+
 #include <future>
 #include <optional>
 #include <queue>
 #include <thread>
 
+#include "debugtrace.hpp"
+#include "lasp_daqdata.h"
+#include "lasp_thread.h"
+
 using namespace std::literals::chrono_literals;
 using lck = std::scoped_lock<std::mutex>;
 using rte = std::runtime_error;
@@ -20,26 +22,26 @@ class SafeQueue {
   std::mutex _mtx;
   std::atomic<uint32_t> _contents{0};
 
-public:
+ public:
   void push(const DaqData &d) {
     DEBUGTRACE_ENTER;
     lck lock(_mtx);
     _queue.push(d);
     _contents++;
-	assert(_contents == _queue.size());
+    assert(_contents == _queue.size());
   }
   DaqData pop() {
     DEBUGTRACE_ENTER;
     if (empty()) {
       throw rte("BUG: Pop on empty queue");
-    } 
+    }
     lck lock(_mtx);
 
     /* DaqData d(std::move(_queue.front())); */
     DaqData d(_queue.front());
     _queue.pop();
     _contents--;
-	assert(_contents == _queue.size());
+    assert(_contents == _queue.size());
     return d;
   }
   /**
@@ -52,58 +54,75 @@ public:
 };
 
 ThreadedInDataHandlerBase::ThreadedInDataHandlerBase(SmgrHandle mgr,
-                                             InCallbackType cb,
-                                             InResetType reset)
-    : _indatahandler(
-          mgr,
-          std::bind(&ThreadedInDataHandlerBase::_inCallbackFromInDataHandler, this,
-                    _1),
-          reset),
-      _queue(std::make_unique<SafeQueue>()), inCallback(cb) {
-
+                                                     InCallbackType cb,
+                                                     ResetCallbackType reset)
+    : _queue(std::make_unique<SafeQueue>()),
+      inCallback(cb),
+      resetCallback(reset),
+      _smgr(mgr) {
   DEBUGTRACE_ENTER;
-
 }
 void ThreadedInDataHandlerBase::startThread() {
   DEBUGTRACE_ENTER;
-  _thread_can_safely_run = true;
-  _indatahandler.start();
+  if (_indatahandler) {
+    throw rte("BUG: ThreadedIndataHandler already started");
+  }
+  SmgrHandle smgr = _smgr.lock();
+  if (!smgr) {
+    cerr << "Stream manager destructed" << endl;
+    return;
+  }
+  _indatahandler = std::make_unique<InDataHandler>(
+      smgr,
+      std::bind(&ThreadedInDataHandlerBase::_inCallbackFromInDataHandler, this,
+                _1),
+      resetCallback);
+
+  _thread_allowed_to_run = true;
+  _indatahandler->start();
 }
 
 void ThreadedInDataHandlerBase::_inCallbackFromInDataHandler(
     const DaqData &daqdata) {
   DEBUGTRACE_ENTER;
-  std::scoped_lock lck(_mtx);
 
   // Early return in case object is under DESTRUCTION
-  if (!_thread_can_safely_run)
-    return;
+  if (!_thread_allowed_to_run) return;
 
   _queue->push(daqdata);
   if (!_thread_running) {
     DEBUGTRACE_PRINT("Pushing new thread in pool");
+    _thread_running = true;
     _pool.push_task(&ThreadedInDataHandlerBase::threadFcn, this);
   }
 }
 
 void ThreadedInDataHandlerBase::stopThread() {
   DEBUGTRACE_ENTER;
-  // Make sure inCallback is no longer called
-  _thread_can_safely_run = false;
-  _indatahandler.stop();
+  if (!_indatahandler) {
+    throw rte("BUG: ThreadedIndataHandler not running");
+  }
 
-  std::scoped_lock lck(_mtx);
+  // Stop the existing thread
+  _thread_allowed_to_run = false;
 
+  // Make sure no new data arrives
+  _indatahandler->stop();
+  _indatahandler.reset();
+
+  DEBUGTRACE_PRINT("Indatahandler stopped. Waiting for thread to finish...");
   // Then wait in steps for the thread to stop running.
   while (_thread_running) {
     std::this_thread::sleep_for(10us);
   }
+  DEBUGTRACE_PRINT("Thread stopped");
+  // Kill the handler
+  DEBUGTRACE_PRINT("Handler resetted");
 }
 
 ThreadedInDataHandlerBase::~ThreadedInDataHandlerBase() {
-
   DEBUGTRACE_ENTER;
-  if (_thread_can_safely_run) {
+  if (_thread_allowed_to_run) {
     stopThread();
     cerr << "*** BUG: InDataHandlers have not been all stopped, while "
             "StreamMgr destructor is called. This is a misuse BUG."
@@ -113,12 +132,9 @@ ThreadedInDataHandlerBase::~ThreadedInDataHandlerBase() {
 }
 
 void ThreadedInDataHandlerBase::threadFcn() {
-
   DEBUGTRACE_ENTER;
-  _thread_running = true;
-
-  while (!_queue->empty() && _thread_can_safely_run) {
 
+  while (!_queue->empty() && _thread_allowed_to_run) {
     // Call inCallback_threaded
     inCallback(_queue->pop());
   }

cpp_src/dsp/lasp_threadedindatahandler.h

@@ -29,21 +29,27 @@ class ThreadedInDataHandlerBase {
    * @brief The queue used to push elements to the handling thread.
    */
 
-  InDataHandler _indatahandler;
   std::unique_ptr<SafeQueue> _queue;
 
-  mutable std::recursive_mutex _mtx;
-
-  std::atomic<bool> _thread_running{false};
-  std::atomic<bool> _thread_can_safely_run{false};
-
-  GlobalThreadPool _pool;
-
   /**
    * @brief Function pointer that is called when new DaqData arrives.
    */
   const InCallbackType inCallback;
 
+  /**
+   * @brief Function pointer that is called when reset() is called.
+   */
+  const ResetCallbackType resetCallback;
+
+  std::weak_ptr<StreamMgr> _smgr;
+  
+  std::unique_ptr<InDataHandler> _indatahandler;
+
+  std::atomic<bool> _thread_running{false};
+  std::atomic<bool> _thread_allowed_to_run{false};
+
+  GlobalThreadPool _pool;
+
   void threadFcn();
 
 
@@ -58,7 +64,7 @@ class ThreadedInDataHandlerBase {
   void _inCallbackFromInDataHandler(const DaqData &daqdata);
 
   public:
-  ThreadedInDataHandlerBase(SmgrHandle mgr, InCallbackType cb, InResetType reset);
+  ThreadedInDataHandlerBase(SmgrHandle mgr, InCallbackType cb, ResetCallbackType reset);
   ~ThreadedInDataHandlerBase();
   /**
    * @brief This method should be called from the derived class' constructor,

cpp_src/dsp/lasp_window.cpp

@@ -11,7 +11,7 @@ using std::cerr;
 using std::endl;
 
 // Safe some typing. Linspace form 0 up to (and NOT including N).
-#define lin0N arma::linspace(0, N - 1, N)
+#define lin0N arma::linspace<vd>(0, N - 1, N)
 
 vd Window::hann(const us N) {
   return arma::pow(arma::sin((arma::datum::pi/N) * lin0N), 2);
@@ -24,8 +24,8 @@ vd Window::blackman(const us N) {
   d a0 = 7938. / 18608.;
   d a1 = 9240. / 18608.;
   d a2 = 1430. / 18608.;
-  return a0 - a1 * d_cos((2 * number_pi/N) * lin0N) +
-         a2 * d_cos((4 * number_pi / N)* lin0N );
+  return a0 - a1 * arma::cos((2 * number_pi/N) * lin0N) +
+         a2 * arma::cos((4 * number_pi / N)* lin0N );
 }
 
 vd Window::rectangular(const us N) { return arma::ones(N); }

cpp_src/lasp_cpp.cpp

@@ -44,6 +44,12 @@ void init_siggen(py::module &m);
 
 PYBIND11_MODULE(lasp_cpp, m) {
 
+#if LASP_DOUBLE_PRECISION == 1
+  m.attr("LASP_DOUBLE_PRECISION") = true;
+#else
+  m.attr("LASP_DOUBLE_PRECISION") = false;
+#endif
+
   init_dsp(m);
   init_deviceinfo(m);
   init_daqconfiguration(m);
@@ -51,6 +57,5 @@ PYBIND11_MODULE(lasp_cpp, m) {
   init_streammgr(m);
   init_datahandler(m);
   init_siggen(m);
-
 }
-/** @} */
+/** @} */

cpp_src/pybind11/lasp_dsp_pybind.cpp

@@ -1,13 +1,16 @@
+#include <pybind11/pybind11.h>
+
+#include <iostream>
+
 #include "arma_npy.h"
 #include "lasp_avpowerspectra.h"
 #include "lasp_biquadbank.h"
 #include "lasp_fft.h"
 #include "lasp_filter.h"
+#include "lasp_freqsmooth.h"
 #include "lasp_slm.h"
 #include "lasp_streammgr.h"
 #include "lasp_window.h"
-#include <iostream>
-#include <pybind11/pybind11.h>
 
 using std::cerr;
 using std::endl;
@@ -27,7 +30,6 @@ using rte = std::runtime_error;
  */
 
 void init_dsp(py::module &m) {
-
   py::class_<Fft> fft(m, "Fft");
   fft.def(py::init<us>());
   fft.def("fft", [](Fft &f, dpyarray dat) {
@@ -114,9 +116,10 @@ void init_dsp(py::module &m) {
 
   aps.def("compute", [](AvPowerSpectra &aps, dpyarray timedata) {
     std::optional<ccube> res;
+    dmat timedata_mat = NpyToMat<d, false>(timedata);
     {
       py::gil_scoped_release release;
-      res = aps.compute(NpyToMat<d, false>(timedata));
+      res = aps.compute(timedata_mat);
     }
 
     return CubeToNpy<c>(res.value_or(ccube(0, 0, 0)));
@@ -151,5 +154,12 @@ void init_dsp(py::module &m) {
   slm.def("Lmax", [](const SLM &slm) { return ColToNpy<d>(slm.Lmax()); });
   slm.def("Lpeak", [](const SLM &slm) { return ColToNpy<d>(slm.Lpeak()); });
   slm.def_static("suggestedDownSamplingFac", &SLM::suggestedDownSamplingFac);
+
+  // Frequency smoother
+  m.def("freqSmooth", [](dpyarray freq, dpyarray X, unsigned w) {
+    vd freqa = NpyToCol<d, false>(freq);
+    vd Xa = NpyToCol<d, false>(X);
+    return ColToNpy(freqSmooth(freqa, Xa, w));
+  });
 }
 /** @} */

cpp_src/pybind11/lasp_pyindatahandler.cpp

@@ -1,4 +1,10 @@
-/* #define DEBUGTRACE_ENABLED */
+// #define DEBUGTRACE_ENABLED
+#include <pybind11/pybind11.h>
+#include <pybind11/pytypes.h>
+
+#include <armadillo>
+#include <atomic>
+
 #include "arma_npy.h"
 #include "debugtrace.hpp"
 #include "lasp_clip.h"
@@ -9,14 +15,12 @@
 #include "lasp_rtsignalviewer.h"
 #include "lasp_streammgr.h"
 #include "lasp_threadedindatahandler.h"
-#include <armadillo>
-#include <atomic>
-#include <chrono>
-#include <pybind11/pybind11.h>
 
 using namespace std::literals::chrono_literals;
 using std::cerr;
 using std::endl;
+using rte = std::runtime_error;
+using Lck = std::scoped_lock<std::recursive_mutex>;
 
 namespace py = pybind11;
 
@@ -48,17 +52,17 @@ py::array_t<T> getPyArrayNoCpy(const DaqData &d) {
    */
 
   return py::array_t<T>(
-      py::array::ShapeContainer({d.nframes, d.nchannels}), // Shape
+      py::array::ShapeContainer({d.nframes, d.nchannels}),  // Shape
 
-      py::array::StridesContainer( // Strides
+      py::array::StridesContainer(  // Strides
           {sizeof(T),
-           sizeof(T) * d.nframes}), // Strides (in bytes) for each index
+           sizeof(T) * d.nframes}),  // Strides (in bytes) for each index
 
       reinterpret_cast<T *>(
-          const_cast<DaqData &>(d).raw_ptr()), // Pointer to buffer
+          const_cast<DaqData &>(d).raw_ptr()),  // Pointer to buffer
 
-      dummyDataOwner // As stated above, now Numpy does not take ownership of
-                     // the data pointer.
+      dummyDataOwner  // As stated above, now Numpy does not take ownership of
+                      // the data pointer.
   );
 }
 
@@ -81,17 +85,17 @@ py::array_t<d> dmat_to_ndarray(const DaqData &d) {
    */
 
   return py::array_t<T>(
-      py::array::ShapeContainer({d.nframes, d.nchannels}), // Shape
+      py::array::ShapeContainer({d.nframes, d.nchannels}),  // Shape
 
-      py::array::StridesContainer( // Strides
+      py::array::StridesContainer(  // Strides
           {sizeof(T),
-           sizeof(T) * d.nframes}), // Strides (in bytes) for each index
+           sizeof(T) * d.nframes}),  // Strides (in bytes) for each index
 
       reinterpret_cast<T *>(
-          const_cast<DaqData &>(d).raw_ptr()), // Pointer to buffer
+          const_cast<DaqData &>(d).raw_ptr()),  // Pointer to buffer
 
-      dummyDataOwner // As stated above, now Numpy does not take ownership of
-                     // the data pointer.
+      dummyDataOwner  // As stated above, now Numpy does not take ownership of
+                      // the data pointer.
   );
 }
 
@@ -104,9 +108,11 @@ class PyIndataHandler : public ThreadedInDataHandler<PyIndataHandler> {
   /**
    * @brief The callback functions that is called.
    */
-  py::function cb, reset_callback;
+  py::object _cb, _reset_callback;
+  std::atomic<bool> _done{false};
+  std::recursive_mutex _mtx;
 
-public:
+ public:
   /**
    * @brief Initialize PyIndataHandler
    *
@@ -117,19 +123,26 @@ public:
    * is called, when a stream stops, this pointer / handle will dangle.
    */
   PyIndataHandler(SmgrHandle mgr, py::function cb, py::function reset_callback)
-      : ThreadedInDataHandler(mgr), cb(cb), reset_callback(reset_callback) {
-
+      : ThreadedInDataHandler(mgr),
+        _cb(py::weakref(cb)),
+        _reset_callback(py::weakref(reset_callback)) {
     DEBUGTRACE_ENTER;
+    // cerr << "Thread ID: " << std::this_thread::get_id() << endl;
     /// Start should be called externally, as at constructor time no virtual
     /// functions should be called.
-    py::gil_scoped_release release;
+    if (_cb().is_none() || _reset_callback().is_none()) {
+      throw rte("cb or reset_callback is none!");
+    }
     startThread();
   }
   ~PyIndataHandler() {
     DEBUGTRACE_ENTER;
+    // cerr << "Thread ID: " << std::this_thread::get_id() << endl;
     /// Callback cannot be called, which results in a deadlock on the GIL
     /// without this release.
     py::gil_scoped_release release;
+    DEBUGTRACE_PRINT("Gil released");
+    _done = true;
     stopThread();
   }
   /**
@@ -137,81 +150,123 @@ public:
    *
    * @param daq Daq device, or nullptr in case no input stream is running.
    */
-  void reset(const Daq *daq) {
+  void reset(const Daq *daqi) {
     DEBUGTRACE_ENTER;
-    try {
-      py::gil_scoped_acquire acquire;
-      if (daq) {
-        reset_callback(daq);
-      } else {
-        reset_callback(py::none());
+    // cerr << "Thread ID: " << std::this_thread::get_id() << endl;
+    if (_done) return;
+    {
+      try {
+        py::object reset_callback = _reset_callback();
+        if (reset_callback.is_none()) {
+          DEBUGTRACE_PRINT("reset_callback is none, weakref killed");
+          _done = true;
+          return;
+        }
+        if (daqi != nullptr) {
+          assert(reset_callback);
+          reset_callback(daqi);
+        } else {
+          assert(reset_callback);
+          reset_callback(py::none());
+        }
+      } catch (py::error_already_set &e) {
+        cerr << "*************** Error calling reset callback!\n";
+        cerr << e.what() << endl;
+        cerr << "*************** \n";
+        /// Throwing a runtime error here does not work out one way or another.
+        /// Therefore, it is better to dive out and prevent undefined behaviour
+        abort();
+        /* throw std::runtime_error(e.what()); */
+      } catch (std::exception &e) {
+        cerr << "Caught unknown exception in reset callback:" << e.what()
+             << endl;
+        abort();
       }
-    } catch (py::error_already_set &e) {
-      cerr << "*************** Error calling reset callback!\n";
-      cerr << e.what() << endl;
-      cerr << "*************** \n";
-      /// Throwing a runtime error here does not work out one way or another.
-      /// Therefore, it is better to dive out and prevent undefined behaviour
-      abort();
-      /* throw std::runtime_error(e.what()); */
-    } catch (std::exception &e) {
-      cerr << "Caught unknown exception in reset callback:" << e.what() << endl;
-      abort();
-    }
-  }
+    }  // end of GIL scope
+  }    // end of function reset()
 
   /**
    * @brief Calls the Python callback method / function with a Numpy array of
    * stream data.
    */
   void inCallback(const DaqData &d) {
-
-    /* DEBUGTRACE_ENTER; */
+    DEBUGTRACE_ENTER;
+    // cerr << "=== Enter incallback for thread ID: " << std::this_thread::get_id() << endl;
 
     using DataType = DataTypeDescriptor::DataType;
-
-    try {
-      py::gil_scoped_acquire acquire;
-      py::object bool_val;
-      switch (d.dtype) {
-      case (DataType::dtype_int8): {
-        bool_val = cb(getPyArrayNoCpy<int8_t>(d));
-      } break;
-      case (DataType::dtype_int16): {
-        bool_val = cb(getPyArrayNoCpy<int16_t>(d));
-      } break;
-      case (DataType::dtype_int32): {
-        bool_val = cb(getPyArrayNoCpy<int32_t>(d));
-      } break;
-      case (DataType::dtype_fl32): {
-        bool_val = cb(getPyArrayNoCpy<float>(d));
-      } break;
-      case (DataType::dtype_fl64): {
-        bool_val = cb(getPyArrayNoCpy<double>(d));
-      } break;
-      default:
-        throw std::runtime_error("BUG");
-      } // End of switch
-
-      bool res = bool_val.cast<bool>();
-    } catch (py::error_already_set &e) {
-      cerr << "ERROR: Python raised exception from callback function: ";
-      cerr << e.what() << endl;
-      abort();
-    } catch (py::cast_error &e) {
-      cerr << e.what() << endl;
-      cerr << "ERROR: Python callback does not return boolean value." << endl;
-      abort();
-    } catch (std::exception &e) {
-      cerr << "Caught unknown exception in Python callback:" << e.what()
-           << endl;
-      abort();
+    if (_done) {
+      DEBUGTRACE_PRINT("Early stop, done");
+      return;
     }
-  }
+    {
+      DEBUGTRACE_PRINT("================ TRYING TO OBTAIN GIL in inCallback...");
+      py::gil_scoped_acquire acquire;
+      try {
+        py::object py_bool;
+        py::object cb = _cb();
+        if (cb.is_none()) {
+          DEBUGTRACE_PRINT("cb is none, weakref killed");
+          _done = true;
+          return;
+        }
+        switch (d.dtype) {
+          case (DataType::dtype_int8): {
+            py_bool = cb(getPyArrayNoCpy<int8_t>(d));
+          } break;
+          case (DataType::dtype_int16): {
+            py_bool = cb(getPyArrayNoCpy<int16_t>(d));
+          } break;
+          case (DataType::dtype_int32): {
+            py_bool = cb(getPyArrayNoCpy<int32_t>(d));
+          } break;
+          case (DataType::dtype_fl32): {
+            py_bool = cb(getPyArrayNoCpy<float>(d));
+          } break;
+          case (DataType::dtype_fl64): {
+            py_bool = cb(getPyArrayNoCpy<double>(d));
+          } break;
+          default:
+            throw std::runtime_error("BUG");
+        }  // End of switch
+
+        bool res = py_bool.cast<bool>();
+        if (res == false) {
+          DEBUGTRACE_PRINT("Setting callbacks to None")
+          _done = true;
+
+          // By doing this, we remove the references, but in the mean time this
+          // might also trigger removing Python objects. Including itself, as
+          // there is no reference to it anymore. The consequence is that the
+          // current object might be destroyed from this thread. However, if we
+          // do not remove these references and in lasp_record.py finish() is
+          // not called, we end up with not-garbage collected recordings in
+          // memory. This is also not good. How can we force Python to not yet
+          // destroy this object?
+          // cb.reset();
+          // reset_callback.reset();
+        }
+      } catch (py::error_already_set &e) {
+        cerr << "ERROR (BUG): Python raised exception from callback function: ";
+        cerr << e.what() << endl;
+        abort();
+      } catch (py::cast_error &e) {
+        cerr << e.what() << endl;
+        cerr << "ERROR (BUG): Python callback does not return boolean value."
+             << endl;
+        abort();
+      } catch (std::exception &e) {
+        cerr << "Caught unknown exception in Python callback:" << e.what()
+             << endl;
+        abort();
+      }
+
+    }  // End of scope in which the GIL is acquired
+    // cerr << "=== LEAVE incallback for thread ID: " << std::this_thread::get_id() << endl;
+
+  }    // End of function inCallback()
 };
 
 void init_datahandler(py::module &m) {
-
   /// The C++ class is PyIndataHandler, but for Python, it is called
   /// InDataHandler
   py::class_<PyIndataHandler> pyidh(m, "InDataHandler");
@@ -244,29 +299,29 @@ void init_datahandler(py::module &m) {
       cval = clip.getCurrentValue();
     }
 
-    return ColToNpy<arma::uword>(cval); // something goes wrong here
+    return ColToNpy<arma::uword>(cval);  // something goes wrong here
   });
 
   /// Real time Aps
   ///
   py::class_<RtAps> rtaps(m, "RtAps");
-  rtaps.def(py::init<SmgrHandle,               // StreamMgr
-                     Filter *const,            // FreqWeighting filter
-                     const us,                 // Nfft
-                     const Window::WindowType, // Window
-                     const d,                  // Overlap percentage 0<=o<100
+  rtaps.def(py::init<SmgrHandle,                // StreamMgr
+                     Filter *const,             // FreqWeighting filter
+                     const us,                  // Nfft
+                     const Window::WindowType,  // Window
+                     const d,                   // Overlap percentage 0<=o<100
 
-                     const d // Time constant
+                     const d  // Time constant
                      >(),
-            py::arg("streammgr"), // StreamMgr
+            py::arg("streammgr"),  // StreamMgr
             py::arg("preFilter").none(true),
             /// Below list of arguments *SHOULD* be same as for
 
             /// AvPowerSpectra constructor!
-            py::arg("nfft") = 2048,                           //
-            py::arg("windowType") = Window::WindowType::Hann, //
-            py::arg("overlap_percentage") = 50.0,             //
-            py::arg("time_constant") = -1                     //
+            py::arg("nfft") = 2048,                            //
+            py::arg("windowType") = Window::WindowType::Hann,  //
+            py::arg("overlap_percentage") = 50.0,              //
+            py::arg("time_constant") = -1                      //
   );
 
   rtaps.def("getCurrentValue", [](RtAps &rt) {
@@ -281,10 +336,10 @@ void init_datahandler(py::module &m) {
   /// Real time Signal Viewer
   ///
   py::class_<RtSignalViewer> rtsv(m, "RtSignalViewer");
-  rtsv.def(py::init<SmgrHandle, // StreamMgr
-                    const d,    // Time history
-                    const us,   // Resolution
-                    const us    // Channel number
+  rtsv.def(py::init<SmgrHandle,  // StreamMgr
+                    const d,     // Time history
+                    const us,    // Resolution
+                    const us     // Channel number
                     >());
 
   rtsv.def("getCurrentValue", [](RtSignalViewer &rt) {

pyproject.toml

@@ -5,7 +5,7 @@ requires-python = ">=3.10"
 description = "Library for Acoustic Signal Processing"
 license = { "file" = "LICENSE" }
 authors = [{ "name" = "J.A. de Jong", "email" = "j.a.dejong@ascee.nl" }]
-version = "1.4.5"
+version = "1.6.6"
 
 keywords = ["DSP", "DAQ", "Signal processing"]
 
@@ -23,8 +23,7 @@ classifiers = [
 urls = { "Documentation" = "https://lasp.ascee.nl" }
 
 dependencies = [
-    "scipy",
-    "numpy",
+    "scipy>=1.13.1",
     "matplotlib>=3.7.2",
     "appdirs",
     "dataclasses_json",
@@ -60,10 +59,3 @@ install_components = ["python_modules"]
 # This might not work properly on Windows. Comment this out when testing on
 # Windows.
 mode = "symlink"
-
-[tool.commitizen]
-name = "cz_conventional_commits"
-tag_format = "v$version"
-version_scheme = "semver"
-version_provider = "pep621"
-update_changelog_on_bump = true

python_src/lasp/filter/fir_design.py

@@ -11,7 +11,9 @@ __all__ = ['freqResponse', 'bandpass_fir_design', 'lowpass_fir_design',
            'arbitrary_fir_design']
 
 import numpy as np
-from scipy.signal import freqz, hann, firwin2
+from scipy.signal import freqz, firwin2
+from scipy.signal.windows import hann
+from ..lasp_config import empty
 
 
 def freqResponse(fs, freq, coefs_b, coefs_a=1.):
@@ -44,7 +46,7 @@ def bandpass_fir_design(L, fs, fl, fu, window=hann):
     Omg2 = 2*np.pi*fu/fs
     Omg1 = 2*np.pi*fl/fs
 
-    fir = np.empty(L, dtype=float)
+    fir = empty(L, dtype=float)
 
     # First Create ideal band-pass filter
     fir[L//2] = (Omg2-Omg1)/np.pi
@@ -64,7 +66,7 @@ def lowpass_fir_design(L, fs, fc, window=hann):
                       " than upper cut-off"
 
     Omgc = 2*np.pi*fc/fs
-    fir = np.empty(L, dtype=float)
+    fir = empty(L, dtype=float)
 
     # First Create ideal band-pass filter
     fir[L//2] = Omgc/np.pi

python_src/lasp/lasp_config.py

@@ -6,8 +6,14 @@ Author: J.A. de Jong - ASCEE
 Description: LASP configuration
 """
 import numpy as np
+from .lasp_cpp import LASP_DOUBLE_PRECISION
 
-LASP_NUMPY_FLOAT_TYPE = np.float64
+if LASP_DOUBLE_PRECISION:
+    LASP_NUMPY_FLOAT_TYPE = np.float64
+    LASP_NUMPY_COMPLEX_TYPE = np.complex128
+else:
+    LASP_NUMPY_FLOAT_TYPE = np.float32
+    LASP_NUMPY_COMPLEX_TYPE = np.float64
 
 
 def zeros(shape):

python_src/lasp/lasp_imptube.py

@@ -1,173 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""!
-Author: J.A. de Jong - ASCEE
-
-Description: Two-microphone impedance tube methods
-"""
-__all__ = ['TwoMicImpedanceTube']
-# from lrftubes import Air
-from .lasp_measurement import Measurement
-from numpy import pi, sqrt, exp
-import numpy as np
-from scipy.interpolate import UnivariateSpline
-# from lrftubes import PrsDuct
-from functools import lru_cache
-
-class TwoMicImpedanceTube:
-    def __init__(self, mnormal: Measurement,
-                       mswitched: Measurement,
-                 s: float,
-                 d1: float,
-                 d2: float,
-                 fl: float = None,
-                 fu: float = None,
-                 periodic_method=False,
-                 # mat= Air(),
-                 D_imptube = 50e-3,
-                 thermoviscous = True,
-                 **kwargs):
-
-        """
-        Initialize two-microphone impedance tube methods
-
-        Args:
-            mnormal: Measurement in normal configuration
-            mswitched: Measurement in normal configuration
-            s: Microphone distance
-            fl: Lower evaluation frequency
-            fu: Lower evaluation frequency
-
-            kwargs: tuple with extra arguments, of which:
-                N: Period length of periodic excitation *obligatory*
-                chan0: Measurement channel index of mic 0
-                chan1: Measurement channel index of mic 1
-
-        """
-        self.mnormal = mnormal
-        self.mswitched = mswitched
-
-        self.mat = mat
-        self.s = s
-        self.d1 = d1
-        self.d2 = d2
-
-        self.fl = fl
-        if fl is None:
-            ksmin = 0.1*pi
-            kmin = ksmin/s
-            self.fl = kmin*mat.c0/2/pi
-
-        self.fu = fu
-        if fu is None:
-            ksmax = 0.8*pi
-            kmax = ksmax/s
-            self.fu = kmax*mat.c0/2/pi
-
-        self.thermoviscous = thermoviscous
-        self.D_imptube = D_imptube
-
-        self.periodic_method = periodic_method
-        self.channels = [kwargs.pop('chan0', 0), kwargs.pop('chan1', 1)]
-        # Compute calibration correction
-        if periodic_method:
-            self.N = kwargs.pop('N')
-            freq, C1 = mnormal.periodicCPS(self.N, channels=self.channels)
-            freq, C2 = mswitched.periodicCPS(self.N, channels=self.channels)
-        else:
-            self.nfft = kwargs.pop('nfft', 16000)
-            freq, C1 = mnormal.CPS(nfft=self.nfft, channels=self.channels)
-            freq, C2 = mswitched.CPS(nfft=self.nfft, channels=self.channels)
-
-        # Store this, as it is often used for just a single sample.
-        self.C1 = C1 
-        self.freq = freq
-
-        self.il = np.where(self.freq<= self.fl)[0][-1]
-        self.ul = np.where(self.freq > self.fu)[0][0]
-
-        # Calibration correction factor
-        # self.K = 0*self.freq + 1.0
-        K = sqrt(C2[:,0,1]*C1[:,0,0]/(C2[:,1,1]*C1[:,1,0]))
-        # self.K = UnivariateSpline(self.freq, K.real)(self.freq) +\
-        #         1j*UnivariateSpline(self.freq, K.imag)(self.freq)
-        self.K = K
-
-    def cut_to_limits(self, ar):
-        return ar[self.il:self.ul]
-
-    def getFreq(self):
-        """
-        Array of frequencies, cut to limits of validity
-        """
-        return self.cut_to_limits(self.freq)
-    
-    @lru_cache
-    def G_AB(self, meas):
-        if meas is self.mnormal:
-            C = self.C1
-            freq = self.freq
-        else:
-            if self.periodic_method:
-                freq, C = meas.periodicCPS(self.N, self.channels)
-            else:
-                freq, C = meas.CPS(nfft=self.nfft, channels=self.channels)
-
-        # Microphone transfer function
-        G_AB = self.K*C[:,1,0]/C[:,0,0]
-        return self.getFreq(), self.cut_to_limits(G_AB)
-    
-    def k(self, freq):
-        """
-        Wave number, or thermoviscous wave number
-        """
-        if self.thermoviscous:
-            D = self.D_imptube
-            S = pi/4*D**2
-            d = PrsDuct(0, S=S, rh=D/4, cs='circ')
-            d.mat = self.mat
-            omg = 2*pi*freq
-            G, Z = d.GammaZc(omg)
-            return G
-        else:
-            return 2*pi*freq/self.mat.c0
-
-
-    def R(self, meas):
-        freq, G_AB = self.G_AB(meas)
-        s = self.s
-        k = self.k(freq)
-        d1 = self.d1
-        RpA = (G_AB - exp(-1j*k*s))/(exp(1j*k*s)-G_AB)
-
-        R = RpA*exp(2*1j*k*(s+d1))
-
-        return freq, R
-
-    def alpha(self, meas):
-        """
-        Acoustic absorption coefficient
-        """
-        freq, R = self.R(meas)
-        return freq, 1 - np.abs(R)**2
-    
-    def z(self, meas):
-        """
-        Acoustic impedance at the position of the sample, in front of the sample
-        """
-        freq, R = self.R(meas)
-        return freq, self.mat.z0*(1+R)/(1-R)
-
-    def zs(self, meas):
-        """
-        Sample impedance jump, assuming a cavity behind the sample with
-        thickness d2
-        """
-        freq, R = self.R(meas)
-        z0 = self.mat.z0
-        k = 2*pi*freq/self.mat.c0
-        d2 = self.d2
-
-        zs = 2*z0*(1-R*exp(2*1j*k*d2))/((R-1)*(exp(2*d2*k*1j)-1))
-        return freq, zs
-

python_src/lasp/lasp_measurement.py

@@ -1,6 +1,24 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 from __future__ import annotations
+from contextlib import contextmanager
+from weakref import WeakValueDictionary
+import h5py as h5
+import uuid
+import pathlib
+import glob
+import itertools
+import numpy as np
+from enum import Enum, unique
+from .lasp_config import LASP_NUMPY_FLOAT_TYPE
+from scipy.io import wavfile
+import os, time, wave, logging
+from .lasp_common import SIQtys, Qty, getFreq
+from .lasp_version import LASP_VERSION_MAJOR, LASP_VERSION_MINOR
+from .lasp_cpp import Window, DaqChannel, AvPowerSpectra
+from typing import List
+from functools import lru_cache
+from .lasp_config import ones
 """!
 Author: J.A. de Jong - ASCEE
 
@@ -46,28 +64,11 @@ The video dataset can possibly be not present in the data.
 
 """
 
-__all__ = ["Measurement", "scaleBlockSens", "MeasurementType"]
-from contextlib import contextmanager
-from weakref import WeakValueDictionary
-import h5py as h5
-import uuid
-import pathlib
-import glob
-import itertools
-import numpy as np
-from enum import Enum, unique
-from .lasp_config import LASP_NUMPY_FLOAT_TYPE
-from scipy.io import wavfile
-import os, time, wave, logging
-from .lasp_common import SIQtys, Qty, getFreq
-from .lasp_version import LASP_VERSION_MAJOR, LASP_VERSION_MINOR
-from .lasp_cpp import Window, DaqChannel, AvPowerSpectra
-from typing import List
-from functools import lru_cache
 
+__all__ = ["Measurement", "scaleBlockSens", "MeasurementType"]
 # Measurement file extension
-MEXT = 'h5'
-DOTMEXT = f'.{MEXT}'
+MEXT = "h5"
+DOTMEXT = f".{MEXT}"
 
 
 @unique
@@ -223,7 +224,7 @@ class IterData(IterRawData):
 
     def __init__(self, fa, channels, sensitivity, **kwargs):
         super().__init__(fa, channels, **kwargs)
-        self.sens = np.asarray(sensitivity)[self.channels]
+        self.sens = np.asarray(sensitivity, dtype=LASP_NUMPY_FLOAT_TYPE)[self.channels]
         assert self.sens.ndim == 1
 
     def __next__(self):
@@ -235,7 +236,8 @@ class Measurement:
     """Provides access to measurement data stored in the h5 measurement file
     format."""
 
-    # Store a dict of open measurements, with uuid string as a key. We store them as a weak ref.
+    # Store a dict of open measurements, with uuid string as a key. We store
+    # them as a weak ref.
     uuid_s = WeakValueDictionary()
 
     def __init__(self, fn):
@@ -271,36 +273,36 @@ class Measurement:
                 self.version_major = f.attrs["LASP_VERSION_MAJOR"]
                 self.version_minor = f.attrs["LASP_VERSION_MINOR"]
             except KeyError:
-                # No version information stored 
+                # No version information stored
                 self.version_major = 0
                 self.version_minor = 1
 
             try:
                 # Try to catch UUID (Universally Unique IDentifier)
-                self._UUID = f.attrs['UUID']
+                self._UUID = f.attrs["UUID"]
                 # Flag indicating we have to add a new UUID
                 create_new_uuid = False
             except KeyError:
                 create_new_uuid = True
 
             try:
-                # UUID of the reference measurement. Should be stored as 
+                # UUID of the reference measurement. Should be stored as
                 # a lists of tuples, where each tuple is a combination of (<MeasurementType.value>, <uuid_string>, <last_filename>).
 
                 # The last filename is a filename that *probably* is the reference measurement with
                 # given UUID. If it is not, we will search for it in the same directory as `this` measurement.
                 # If we cannot find it there, we will give up, and remove the field corresponding to this reference measurement type.
-                refMeas_list = f.attrs['refMeas']
+                refMeas_list = f.attrs["refMeas"]
 
                 # Build a tuple string from it
                 self._refMeas = {}
-                for (key, val, name) in refMeas_list:
+                for key, val, name in refMeas_list:
                     self._refMeas[MeasurementType(int(key))] = (val, name)
             except KeyError:
                 self._refMeas = {}
 
             try:
-                self._type_int = f.attrs['type_int']
+                self._type_int = f.attrs["type_int"]
             except KeyError:
                 self._type_int = 0
 
@@ -329,10 +331,10 @@ class Measurement:
             try:
                 sens = f.attrs["sensitivity"]
                 self._sens = (
-                    sens * np.ones(self.nchannels) if isinstance(sens, float) else sens
+                    sens * ones(self.nchannels) if isinstance(sens, float) else sens
                 )
             except KeyError:
-                self._sens = np.ones(self.nchannels)
+                self._sens = ones(self.nchannels)
 
             # The time is cached AND ALWAYS ASSUMED TO BE AN IMMUTABLE OBJECT.
             # It is also cached. Changing the measurement timestamp should not
@@ -367,7 +369,9 @@ class Measurement:
             self.genNewUUID()
         else:
             if self.UUID in Measurement.uuid_s.keys():
-                raise RuntimeError(f"Measurement '{self.name}' is already opened. Cannot open measurement twice. Note: this error can happen when measurements are manually copied.")
+                raise RuntimeError(
+                    f"Measurement '{self.name}' is already opened. Cannot open measurement twice. Note: this error can happen when measurements are manually copied."
+                )
 
         # Store weak reference to 'self' in list of UUID's. They are removed when no file is open anymore
         Measurement.uuid_s[self._UUID] = self
@@ -379,11 +383,11 @@ class Measurement:
         Args:
             newname: New name, with or without extension
         """
-        _ , ext = os.path.splitext(newname)
+        _, ext = os.path.splitext(newname)
         # Add proper extension if new name is given without extension.
         if ext != DOTMEXT:
             newname = newname + DOTMEXT
-            
+
         # Folder, Base filename + extension
         folder, _ = os.path.split(self.fn)
 
@@ -396,7 +400,7 @@ class Measurement:
         """
         Create new UUID for measurement and store in file.
         """
-        self.setAttribute('UUID', str(uuid.uuid1()))
+        self.setAttribute("UUID", str(uuid.uuid1()))
 
     @property
     def UUID(self):
@@ -404,59 +408,73 @@ class Measurement:
         Universally unique identifier
         """
         return self._UUID
-    
+
     def getRefMeas(self, mtype: MeasurementType):
         """
-        Return corresponding reference measurement, if configured and can be found. If the reference 
+        Return corresponding reference measurement, if configured and can be found. If the reference
         measurement is currently not open, it tries to open it by traversing other measurement
         files in the current directory. Throws a runtime error in case the reference measurement cannot be found.
         Throws a ValueError when the reference measurement is not configured.
         """
-        
+
         # See if we can find the UUID for the required measurement type
         try:
             required_uuid, possible_name = self._refMeas[mtype]
         except KeyError:
             raise ValueError(f"No reference measurement configured for '{self.name}'")
-            
+
         m = None
         # Try to find it in the dictionary of of open measurements
         if required_uuid in Measurement.uuid_s.keys():
             m = Measurement.uuid_s[required_uuid]
-            logging.info(f'Returned reference measurement {m.name} from list of open measurements')
-        
+            logging.debug(
+                f"Returned reference measurement {m.name} from list of open measurements"
+            )
+
         # Not found in list of openend measurements. See if we can open it using its last stored file name we know of
         if m is None:
             try:
                 m = Measurement(possible_name)
                 if m.UUID == required_uuid:
-                    logging.info(f'Opened reference measurement {m.name} by name')
+                    logging.info(f"Opened reference measurement {m.name} by name")
             except Exception as e:
-                logging.error(f'Could not find reference measurement using file name: {possible_name}')
+                logging.error(
+                    f"Could not find reference measurement using file name: {possible_name}"
+                )
 
         # Last resort, see if we can find the right measurement in the same folder
         if m is None:
             try:
                 folder, _ = os.path.split(self.fn)
-                m = Measurement.fromFolderWithUUID(required_uuid, folder, skip=[self.name])
-                logging.info('Found reference measurement in folder with correct UUID. Updating name of reference measurement')
+                m = Measurement.fromFolderWithUUID(
+                    required_uuid, folder, skip=[self.name]
+                )
+                logging.info(
+                    "Found reference measurement in folder with correct UUID. Updating name of reference measurement"
+                )
                 # Update the measurement file name in the list, such that next time it
                 # can be opened just by its name.
                 self.setRefMeas(m)
             except:
-                logging.error("Could not find the reference measurement. Is it deleted?")
+                logging.error(
+                    "Could not find the reference measurement. Is it deleted?"
+                )
 
         # Well, we found it. Now make sure the reference measurement actually has the right type (User could have marked it as a NotSpecific for example in the mean time).
         if m is not None:
             if m.measurementType() != mtype:
                 m.removeRefMeas(mtype)
-                raise RuntimeError(f"Reference measurement for {self.name} is not a proper reference (anymore).")
-            
+                raise RuntimeError(
+                    f"Reference measurement for {self.name} is not a proper reference (anymore)."
+                )
+
             # Whow, we passed all security checks, here we go!
             return m
         else:
             # Nope, not there.
-            raise RuntimeError(f"Could not find the reference measurement for '{self.name}'. Is it deleted?")
+            raise RuntimeError(
+                f"Could not find the reference measurement for '{self.name}'. Is it deleted?"
+            )
 
     def removeRefMeas(self, mtype: MeasurementType):
         """
@@ -477,9 +495,12 @@ class Measurement:
         with self.file("r+") as f:
             # Update attribute in file. Stored as a flat list of string tuples:
             # [(ref_value1, uuid_1, name_1), (ref_value2, uuid_2, name_2), ...]
-            reflist = list((str(key.value), val1, val2) for key, (val1, val2) in self._refMeas.items())
+            reflist = list(
+                (str(key.value), val1, val2)
+                for key, (val1, val2) in self._refMeas.items()
+            )
             # print(reflist)
-            f.attrs['refMeas'] = reflist
+            f.attrs["refMeas"] = reflist
 
     def setRefMeas(self, m: Measurement):
         """
@@ -489,32 +510,36 @@ class Measurement:
         """
         mtype = m.measurementType()
         if mtype == MeasurementType.NotSpecific:
-            raise ValueError('Measurement to be set as reference is not a reference measurement')
-    
+            raise ValueError(
+                "Measurement to be set as reference is not a reference measurement"
+            )
+
         self._refMeas[mtype] = (m.UUID, m.name)
         self.__storeReafMeas()
 
     @staticmethod
-    def fromFolderWithUUID(uuid_str: str, folder: str='', skip=[]):
+    def fromFolderWithUUID(uuid_str: str, folder: str = "", skip=[]):
         """
         Returns Measurement object from a given UUID string. It first tries to find whether there
-        is an uuid in the static list of weak references. If not, it will try to open files in 
+        is an uuid in the static list of weak references. If not, it will try to open files in
         the current file path.
         """
-        for fn in glob.glob(str(pathlib.Path(folder)) + f'/*{DOTMEXT}'):
+        for fn in glob.glob(str(pathlib.Path(folder)) + f"/*{DOTMEXT}"):
             # Do not try to open this file in case it is in the 'skip' list.
             if len(list(filter(lambda a: a in fn, skip))) > 0:
                 continue
-                
+
             try:
                 m = Measurement(fn)
-                if m.UUID == uuid_str: 
+                if m.UUID == uuid_str:
                     # Update 'last_fn' attribute in dict of stored reference measurements
                     return m
             except Exception as e:
-                logging.error(f'Possible measurement file {fn} returned error {e} when opening.')
+                logging.error(
+                    f"Possible measurement file {fn} returned error {e} when opening."
+                )
 
-        raise RuntimeError(f'Measurement with UUID {uuid_str} could not be found.')
+        raise RuntimeError(f"Measurement with UUID {uuid_str} could not be found.")
 
     def setAttribute(self, attrname: str, value):
         """
@@ -534,17 +559,17 @@ class Measurement:
         """
         Returns True when a measurement is flagged as being of a certaint "MeasurementType"
         """
-        if (type_.value & self._type_int):
+        if type_.value & self._type_int:
             return True
         elif type_.value == self._type_int == 0:
             return True
         return False
-    
+
     def setType(self, type_: MeasurementType):
         """
         Set the measurement type to given type
         """
-        self.setAttribute('type_int', type_.value)
+        self.setAttribute("type_int", type_.value)
 
     def measurementType(self):
         """
@@ -589,7 +614,7 @@ class Measurement:
     @channelConfig.setter
     def channelConfig(self, chcfg: List[DaqChannel]):
         """
-        Set new channel configuration from list of DaqChannel objects. 
+        Set new channel configuration from list of DaqChannel objects.
 
         Use cases:
             - Update channel types, sensitivities etc.
@@ -885,7 +910,7 @@ class Measurement:
         """
         if isinstance(sens, float):
             # Put all sensitivities equal
-            sens = sens * np.ones(self.nchannels)
+            sens = sens * ones(self.nchannels)
         elif isinstance(sens, list):
             sens = np.asarray(sens)
 
@@ -998,9 +1023,9 @@ class Measurement:
         happen that a Measurement object is created twice for the same backing file, which we do not allow.
         """
         # See if the base part of the filename is referring to a file that is already open
-        with h5.File(fn, 'r') as f:
+        with h5.File(fn, "r") as f:
             try:
-                theuuid = f.attrs['UUID']
+                theuuid = f.attrs["UUID"]
             except KeyError:
                 # No UUID stored in measurement. This is an old measurement that did not have UUID's
                 # We create a new UUID here such that the file is opened from the filesystem
@@ -1009,7 +1034,7 @@ class Measurement:
 
         if theuuid in Measurement.uuid_s.keys():
             return Measurement.uuid_s[theuuid]
-        
+
         return Measurement(fn)
 
     @staticmethod
@@ -1090,8 +1115,8 @@ class Measurement:
         sensitivity,
         mfn,
         timestamp=None,
-        qtys: List[SIQtys]=None,
-        channelNames: List[str]=None,
+        qtys: List[SIQtys] = None,
+        channelNames: List[str] = None,
         force=False,
     ) -> Measurement:
         """
@@ -1161,7 +1186,7 @@ class Measurement:
             hf.attrs["nchannels"] = nchannels
 
             # Add physical quantity indices
-            hf.attrs['qtys_enum_idx'] = [qty.toInt() for qty in qtys]
+            hf.attrs["qtys_enum_idx"] = [qty.toInt() for qty in qtys]
 
             # Add channel names in case given
             if channelNames is not None:
@@ -1199,7 +1224,7 @@ class Measurement:
             nchannels = 1
             nframes = len(data)
             data = data[:, np.newaxis]
-        sensitivity = np.ones(nchannels)
+        sensitivity = ones(nchannels)
 
         with h5.File(newfn, "w") as hf:
             hf.attrs["samplerate"] = samplerate
@@ -1217,4 +1242,3 @@ class Measurement:
             ad[0] = data
 
         return Measurement(newfn)
-

python_src/lasp/lasp_measurementset.py

@@ -3,7 +3,8 @@ Provides class MeasurementSet, a class used to perform checks and adjustments
 on a group of measurements at the same time.
 
 """
-__all__ = ['MeasurementSet']
+
+__all__ = ["MeasurementSet"]
 from .lasp_measurement import Measurement, MeasurementType
 from typing import List
 import time
@@ -13,28 +14,35 @@ class MeasurementSet(list):
     """
     Group of measurements that have some correspondence to one another. Class
     is used to operate on multiple measurements at once.
-
     """
 
-    def __init__(self, mlist: List[Measurement]=[]):
+    def __init__(self, mlist: List[Measurement] = []):
         """
         Initialize a measurement set
 
-        Args:
+        Arg:
             mlist: Measurement list
-
         """
-        if any([not isinstance(i, Measurement) for i in mlist]):
-            raise TypeError('Object in list should be of Measurement type')
 
+        if any([not isinstance(i, Measurement) for i in mlist]):
+            raise TypeError("Object in list should be of Measurement type")
+
+        # Sort by time stamp, otherwise the order is random
+        mlist.sort(key=lambda x: x.time, reverse=True)
         super().__init__(mlist)
 
     def getNewestReferenceMeasurement(self, mtype: MeasurementType):
-        """Return the newest (in time) measurement in the current list of a certain type. Returns None in case no measurement could be found.
-
-        Args:
-            mtype (MeasurementType): The type required.
         """
+        Get the NEWEST ref. measurement of a current type, in the current set.
+
+        Arg:
+            mtype (MeasurementType): The type required.
+
+        Return:
+            - The newest (in time) measurement in the current list of a certain type.
+            - None, in case no measurement could be found.
+        """
+
         mnewest = None
         for m in self:
             if m.measurementType() == mtype:
@@ -44,15 +52,34 @@ class MeasurementSet(list):
                     if mnewest.time < m.time:
                         mnewest = m
         return mnewest
-    
+
+    def getReferenceMeasurements(self, mtype: MeasurementType):
+        """Get ALL ref. measurements of a certain type, in the current set.
+
+        Arg:
+            mtype (MeasurementType): The type of which to list
+
+        Return:
+            A new measurement set, including all measurements of a certain type
+        """
+
+        return [m for m in self if m.measurementType() == mtype]
+
     def getNewestReferenceMeasurements(self):
-        """Returns a dictionary with newest measurement of each type that is not specific returns None in case no measurement is found."""
+        """
+        Get the NEWEST ref. measurement of all types, in the current set.
+
+        Return:
+            - A dictionary with the newest measurement of each type that is not specific
+            - None, in case no measurement is found.
+        """
+
         newest = {}
         for m in self:
             mtype = m.measurementType()
             if mtype == MeasurementType.NotSpecific:
                 continue
-            if not mtype in newest:
+            if mtype not in newest:
                 newest[mtype] = m
             else:
                 if m.time > newest[mtype].time:
@@ -60,38 +87,50 @@ class MeasurementSet(list):
         return newest
 
     def newestReferenceOlderThan(self, secs):
-        """Returns a dictionary of references with the newest reference, that is still
-        older than `secs` seconds. """
+        """
+        Get a dictionary of reference measurements which are older than a
+        specified threshold. Only one of each type is returned.
+
+        Args:
+            - secs: time threshold, in seconds
+
+        Return:
+            - a dictionary of references with the newest reference, that is still
+            older than `secs` seconds
+        """
         curtime = time.time()
         newest = self.getNewestReferenceMeasurements()
-        newest_older_than =  {}
+        newest_older_than = {}
         for key, m in newest.items():
             if curtime - m.time >= secs:
                 newest_older_than[key] = m
         return newest_older_than
 
-
     def measTimeSame(self):
         """
-        Returns True if all measurements have the same measurement
-        time (recorded time)
-
+        Returns True if all measurements have the same measurement length and
+        sample rate
         """
+
         if len(self) > 0:
-            first = self[0].N
-            return all([first == meas.N for meas in self])
+            firstN = self[0].N              # samples
+            firstFS = self[0].samplerate    # sample rate
+            sameN = all([firstN == meas.N for meas in self])
+            sameFS = all([firstFS == meas.samplerate for meas in self])
+            return sameN and sameFS
         else:
             return False
 
     def measSimilar(self):
         """
         Similar means: channel metadata is the same, and the measurement time
-        is the same. It means that the recorded data is, of course, different.
+        is the same.
 
         Returns:
-            True if measChannelsSame() and measTimeSame() else False
-
+            - True if measChannelsSame() and measTimeSame()
+            - False otherwise
         """
+
         return self.measTimeSame() and self.measChannelsSame()
 
     def measChannelsSame(self):
@@ -101,9 +140,9 @@ class MeasurementSet(list):
         a set of measurements, simultaneously. If the channel data is the same
         (name, sensitivity, ...) it returns True.
         """
+
         if len(self) > 0:
             first = self[0].channelConfig
             return all([first == meas.channelConfig for meas in self])
         else:
             return False
-

python_src/lasp/lasp_octavefilter.py

@@ -77,6 +77,7 @@ class FirFilterBank:
 
         self.fs  = fs
         self.xs = list(range(xmin, xmax + 1))
+        raise RuntimeError('Not working code anymore')
 
         maxdecimation = self.designer.firDecimation(self.xs[0])
         self.decimators = []
@@ -245,7 +246,7 @@ class SosFilterBank:
         for i, x in enumerate(self.xs):
             channel = self.designer.createSOSFilter(x)
             if sos is None:
-                sos = np.empty((channel.size, len(self.xs)))
+                sos = empty((channel.size, len(self.xs)))
             sos[:, i] = channel.flatten() 
 
         self._fb = BiquadBank(sos)

python_src/lasp/lasp_record.py

@@ -11,6 +11,11 @@ from enum import Enum, auto, unique
 from .lasp_cpp import InDataHandler, StreamMgr
 from .lasp_version import LASP_VERSION_MAJOR, LASP_VERSION_MINOR
 import uuid
+import logging
+
+logger = logging.getLogger(__name__)
+# logger.setLevel(logging.DEBUG)
+logger.setLevel(logging.INFO)
 
 
 @dataclasses.dataclass
@@ -60,10 +65,14 @@ class Recording:
             startDelay: Optional delay added before the recording is *actually*
             started in [s].
         """
+        logger.debug("__init__()")
         ext = ".h5"
         if ext not in fn:
             fn += ext
 
+        if os.path.exists(fn):
+            raise RuntimeError("Recording file name already exists / is in use")
+
         self._smgr = streammgr
         self._metadata = None
 
@@ -91,7 +100,7 @@ class Recording:
         self._stop = Atomic(False)
 
         # Mutex, on who is working with the H5py data and the class settings
-        self._rec_mutex = threading.Lock()
+        self._rec_mutex = threading.RLock()
 
         self._progressCallback = progressCallback
 
@@ -100,7 +109,7 @@ class Recording:
             self._h5file = h5py.File(self._fn, "w", "stdio")
             self._h5file.flush()
         except Exception as e:
-            logging.error(f"Error creating measurement file {e}")
+            logger.error(f"Error creating measurement file {e}")
             raise
 
         # This flag is used to delete the file on finish(), and can be used
@@ -115,19 +124,26 @@ class Recording:
         # Audio dataset
         self._ad = None
 
-        logging.debug("Starting record....")
+        logger.debug("Starting record....")
 
+        # In the PyInDataHandler, a weak reference is stored to the python
+        # methods reset and incallback. One way or another, the weak ref is gone
+        # on the callback thread. If we store an "extra" ref to this method over
+        # here, the weak ref stays alive. We do not know whether this is a bug
+        # or a feature, but in any case storing this extra ref to inCallback
+        # solves the problem.
+        self._incalback_cpy = self.inCallback
         self._indataHandler = InDataHandler(
-            streammgr, self.inCallback, self.resetCallback
+            streammgr, self._incalback_cpy, self.resetCallback
         )
 
         if wait:
-            logging.debug("Stop recording with CTRL-C")
+            logger.debug("Stop recording with CTRL-C")
             try:
                 while not self._stop():
                     time.sleep(0.01)
             except KeyboardInterrupt:
-                logging.debug("Keyboard interrupt on record")
+                logger.debug("Keyboard interrupt on record")
             finally:
                 self.finish()
 
@@ -138,6 +154,7 @@ class Recording:
         """
         Function called with initial stream data.
         """
+        logger.debug(f"resetCallback({daq})")
         with self._rec_mutex:
             in_ch = daq.enabledInChannels()
             blocksize = daq.framesPerBlock()
@@ -180,10 +197,11 @@ class Recording:
         When returning False, it will stop the stream.
 
         """
+        logger.debug(f"inCallback()")
         if self._stop():
-            logging.debug("Stop flag set, early return in inCallback")
+            logger.debug("Stop flag set, early return in inCallback")
             # Stop flag is raised. We do not add any data anymore.
-            return True
+            return False
 
         with self._rec_mutex:
 
@@ -191,10 +209,7 @@ class Recording:
 
             match self._recState:
                 case RecordingState.Waiting:
-                    if (
-                        self._allBlocks * self._blocksize / self._fs
-                        > self._startDelay
-                    ):
+                    if self._allBlocks * self._blocksize / self._fs > self._startDelay:
                         self._recState = RecordingState.Recording
 
                 case RecordingState.Recording:
@@ -208,7 +223,10 @@ class Recording:
 
                     recstatus = RecordStatus(curT=self.recordedTime, done=False)
 
-                    if self._requiredRecordingLength is not None and self.recordedTime >= self._requiredRecordingLength:
+                    if (
+                        self._requiredRecordingLength is not None
+                        and self.recordedTime >= self._requiredRecordingLength
+                    ):
                         self._recState = RecordingState.AllDataStored
                         self._stop <<= True
                         recstatus.done = True
@@ -217,19 +235,20 @@ class Recording:
                         self._progressCallback(recstatus)
 
                 case RecordingState.AllDataStored:
-                    pass
+                    return False
 
                 case RecordingState.Finished:
-                    pass
+                    return False
 
             return True
 
     @property
     def recordedTime(self):
         """Return recorded time (not rounded) as float"""
-        if self._ad is None:
-            return 0.0
-        return self._recordedBlocks * self._blocksize / self._fs
+        with self._rec_mutex:
+            if self._ad is None:
+                return 0.0
+            return self._recordedBlocks * self._blocksize / self._fs
 
     def __addFirstFramesToFile(self, adata):
         """
@@ -240,29 +259,29 @@ class Recording:
             adata: Numpy array with data from DAQ
 
         """
+        with self._rec_mutex:
+            # The array data type cannot
+            # datatype = daq.dataType()
+            dtype = np.dtype(adata.dtype)
 
-        # The array data type cannot
-        # datatype = daq.dataType()
-        dtype = np.dtype(adata.dtype)
+            assert self._ad is None
 
-        assert self._ad is None
+            self._ad = self._h5file.create_dataset(
+                "audio",
+                (1, self._blocksize, self._nchannels),
+                dtype=dtype,
+                maxshape=(
+                    None,  # This means, we can add blocks
+                    # indefinitely
+                    self._blocksize,
+                    self._nchannels,
+                ),
+                compression="gzip",
+            )
+            self._ad[0, :, :] = adata
 
-        self._ad = self._h5file.create_dataset(
-            "audio",
-            (1, self._blocksize, self._nchannels),
-            dtype=dtype,
-            maxshape=(
-                None,  # This means, we can add blocks
-                # indefinitely
-                self._blocksize,
-                self._nchannels,
-            ),
-            compression="gzip",
-        )
-        self._ad[0, :, :] = adata
-
-        self._h5file.flush()
-        self._deleteFile = False
+            self._h5file.flush()
+            self._deleteFile = False
 
     def setDelete(self, val: bool):
         """
@@ -279,36 +298,35 @@ class Recording:
         remove the queue from the stream, etc.
 
         """
-        logging.debug("Recording::finish()")
+        logger.debug("Recording::finish()")
 
         self._stop <<= True
 
-
         with self._rec_mutex:
-            if self._recState == RecordingState.Finished:
-                raise RuntimeError('Recording has already finished')
 
-            self._h5file.flush()
+            if self._recState == RecordingState.Finished:
+                raise RuntimeError("Recording has already finished")
+
             # Remove indata handler, which also should remove callback function
             # from StreamMgr. This, however does not have to happen
             # instantaneously. For which we have to implement extra mutex
             # guards in this class
             del self._indataHandler
-            self._indataHandler = None
+
+            self._h5file.flush()
 
             # Remove handle to dataset otherwise the h5 file is not closed
             # properly.
             del self._ad
-            self._ad = None
 
             try:
                 # Close the recording file
                 self._h5file.close()
                 del self._h5file
             except Exception as e:
-                logging.error(f"Error closing file: {e}")
+                logger.error(f"Error closing file: {e}")
 
-            logging.debug("Recording ended")
+            logger.debug("Recording ended")
             if self._deleteFile:
                 self.__deleteFile()
             self._recState = RecordingState.Finished
@@ -320,16 +338,17 @@ class Recording:
         try:
             os.remove(self._fn)
         except Exception as e:
-            logging.error(f"Error deleting file: {self._fn}: {str(e)}")
+            logger.error(f"Error deleting file: {self._fn}: {str(e)}")
 
     def __addTimeDataToFile(self, indata):
         """
         Called by handleQueue() and adds new time data to the storage file.
         """
+        with self._rec_mutex:
 
-        ablockno = self._recordedBlocks
+            ablockno = self._recordedBlocks
 
-        # Add the data to the file, and resize the audio data blocks
-        self._ad.resize(ablockno + 1, axis=0)
-        self._ad[ablockno, :, :] = indata
-        self._h5file.flush()
+            # Add the data to the file, and resize the audio data blocks
+            self._ad.resize(ablockno + 1, axis=0)
+            self._ad[ablockno, :, :] = indata
+            self._h5file.flush()

python_src/lasp/lasp_reverb.py

@@ -6,6 +6,7 @@ Description:
 Reverberation time estimation tool using least squares
 """
 from .lasp_common import getTime
+from .lasp_config import ones
 import numpy as np
 
 
@@ -56,7 +57,7 @@ class ReverbTime:
         x = self._t[istart:istop][:, np.newaxis]
 
         # Solve the least-squares problem, by creating a matrix of
-        A = np.hstack([x, np.ones(x.shape)])
+        A = np.hstack([x, ones(x.shape)])
 
         # print(A.shape)
         # print(points.shape)

python_src/lasp/lasp_slm.py

@@ -5,6 +5,7 @@ Sound level meter implementation
 @author: J.A. de Jong - ASCEE
 """
 from .lasp_cpp import cppSLM
+from .lasp_config import empty
 import numpy as np
 from .lasp_common import (TimeWeighting, FreqWeighting, P_REF)
 from .filter import SPLFilterDesigner
@@ -101,7 +102,7 @@ class SLM:
             assert fbdesigner.fs == fs
             sos_firstx = fbdesigner.createSOSFilter(self.xs[0]).flatten()
             self.nom_txt.append(fbdesigner.nominal_txt(self.xs[0]))
-            sos = np.empty((sos_firstx.size, nfilters), dtype=float, order='C')
+            sos = empty((sos_firstx.size, nfilters), dtype=float, order='C')
             sos[:, 0] = sos_firstx
 
             for i, x in enumerate(self.xs[1:]):

python_src/lasp/lasp_weighcal.py

@@ -47,7 +47,7 @@ class WeighCal:
 
         P = 2048 # Filter length (number of taps)
 
-        self._firs = np.empty((P, self.nchannels))
+        self._firs = empty((P, self.nchannels))
         self._fbs = []
         for chan in range(self.nchannels):
             fir = arbitrary_fir_design(fs, P, freq_design,

python_src/lasp/tools/tools.py

@@ -20,6 +20,8 @@ from enum import Enum, unique
 import copy
 import numpy as np
 from numpy import log2, pi, sin
+from ..lasp_cpp import freqSmooth
+from ..lasp_config import zeros
 
 
 @unique
@@ -152,7 +154,7 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
     return Q
 
 
-def smoothSpectralData(freq, M, sw: SmoothingWidth,
+def smoothSpectralData_old(freq, M, sw: SmoothingWidth,
                        st: SmoothingType = SmoothingType.levels):
     """
     Apply fractional octave smoothing to data in the frequency domain.
@@ -220,6 +222,45 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 
     return Psm
 
+def smoothSpectralData(freq, M, sw: SmoothingWidth,
+                       st: SmoothingType = SmoothingType.levels):
+    """
+    Apply fractional octave smoothing to data in the frequency domain.
+
+    Args:
+        freq: array of frequencies of data points [Hz] - equally spaced
+        M: array of data, either power or dB
+        the smoothing type `st`, the smoothing is applied.
+        sw: smoothing width
+        st: smoothing type = data type of input data
+
+    Returns:
+        freq : array frequencies of data points [Hz]
+        Msm : float smoothed magnitude of data points
+    """
+    # Safety
+
+    if st == SmoothingType.ps:
+        assert np.min(M) >= 0, 'Power spectrum values cannot be negative'
+    if st == SmoothingType.levels and isinstance(M.dtype, complex):
+        raise RuntimeError('Decibel input should be real-valued')
+
+    # Convert to power
+    if st == SmoothingType.levels:
+        P = 10**(M/10)
+    elif st == SmoothingType.ps:
+        P = M
+    else:
+        raise RuntimeError(f"Incorrect SmoothingType: {st}")
+
+    Psm = freqSmooth(freq, P, sw.value[0])
+
+    # Convert to original format
+    if st == SmoothingType.levels:
+        Psm = 10*np.log10(Psm)
+
+    return Psm
+
 
 # %% Test
 if __name__ == "__main__":

scripts/copy_windows_dlls.sh

@@ -1,11 +1,9 @@
 #!/bin/bash
 #
-cmake . -G"Ninja" -DLASP_HAS_ULDAQ=OFF -DPython3_ROOT_DIR=C:\\winpython\\python-3.10.9.amd64
-arch=ucrt64
-#arch=mingw64
-
+arch=mingw64
+# DLL's that are required by lasp_cpp
 files_to_cpy="libfftw3-3.dll libgcc_s_seh-1.dll libgfortran-5.dll libgomp-1.dll libopenblas.dll libquadmath-0.dll libstdc++-6.dll libwinpthread-1.dll"
 
 for fn in ${files_to_cpy}; do
-	cp /c/msys64/${arch}/bin/${fn} src/lasp
+	cp /c/msys64/${arch}/bin/${fn} python_src/lasp
 done

scripts/install_msys2_builddeps.sh

@@ -5,8 +5,8 @@
 if [ -z $CI ]; then
     PACMAN_OPTIONS="--needed --noconfirm"
 fi
-# arch=mingw-w64-x86_64
-arch=mingw-w64-ucrt-x86_64
+arch=mingw-w64-x86_64
+# arch=mingw-w64-ucrt-x86_64
 
 pacman -S ${PACMAN_OPTIONS} make
 

test/test_fft.py

@@ -31,7 +31,7 @@ def test_backward_fft():
     nfft = 2048
     freq = getFreq(nfft, nfft)
     
-    # Sig = np.zeros(nfft//2+1, dtype=complex)
+    # Sig = zeros(nfft//2+1, dtype=complex)
     Sigr = np.random.randn(nfft//2+1)
     Sigi = np.random.randn(nfft//2+1)
     

third_party/portaudio

@@ -1 +1 @@
-Subproject commit daaf637f6f9fce670031221abfd7dfde92e5cce3
+Subproject commit 18a606e1f928852bfc29639d9539ae74d37b5dee