From 70891ceaf4cdca87c1eb653234ec85e5660e6c22 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Sat, 3 Sep 2022 20:59:14 +0200 Subject: [PATCH] Switched to OpenMP for parallellizing for loops. Bugfixes in PowerSpectra() class. Added tests to check Parsevall hold. Improved lots of comments. Added module groups. Use CMake to specify code version number. Device info should be obtained from StreamMgr in Python code. --- CMakeLists.txt | 63 ++++++--- examples/Test SLM.ipynb | 57 ++------ examples/Test SeriesBiquad.ipynb | 25 ++-- examples/Test input.ipynb | 65 +++++---- pyproject.toml | 1 - src/lasp/CMakeLists.txt | 1 - src/lasp/device/lasp_daq.h | 19 ++- src/lasp/device/lasp_daqconfig.h | 46 ++++-- src/lasp/device/lasp_daqdata.h | 5 + src/lasp/device/lasp_deviceinfo.h | 4 + src/lasp/device/lasp_streammgr.cpp | 86 ++++++++++-- src/lasp/device/lasp_streammgr.h | 64 +++++++-- src/lasp/dsp/lasp_avpowerspectra.cpp | 13 +- src/lasp/dsp/lasp_avpowerspectra.h | 4 + src/lasp/dsp/lasp_biquadbank.cpp | 12 +- src/lasp/dsp/lasp_fft.cpp | 140 +++++++++++-------- src/lasp/dsp/lasp_fft.h | 11 +- src/lasp/dsp/lasp_siggen_impl.h | 3 +- src/lasp/dsp/lasp_slm.cpp | 31 +--- src/lasp/dsp/lasp_slm.h | 2 +- src/lasp/dsp/lasp_threadedindatahandler.h | 4 + src/lasp/dsp/lasp_timebuffer.cpp | 6 + src/lasp/dsp/lasp_types.h | 25 ++-- src/lasp/lasp_config.h.in | 4 +- src/lasp/lasp_config.py | 6 +- src/lasp/lasp_slm.py | 4 +- src/lasp/pybind11/lasp_deviceinfo.cpp | 2 - src/lasp/pybind11/lasp_dsp_pybind.cpp | 22 ++- src/lasp/pybind11/lasp_streammgr.cpp | 1 + test/test_biquadbank.py | 40 ++++++ test/test_cppslm.py | 75 ++++++++++ test/test_fft.c | 33 ----- test/test_fft.py | 3 + test/test_math.c | 50 ------- test/test_ps.py | 163 +++++++++++++++++----- test/test_rtaudio.py | 35 ----- 36 files changed, 699 insertions(+), 426 deletions(-) create mode 100644 test/test_biquadbank.py create mode 100644 test/test_cppslm.py delete mode 100644 test/test_fft.c delete mode 100644 test/test_math.c delete mode 100644 test/test_rtaudio.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 74ef5a5..7ec236e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,22 +1,6 @@ cmake_minimum_required (VERSION 3.16) +project(LASP LANGUAGES C CXX VERSION 1.0) -project(LASP LANGUAGES C CXX) -set(LASP_VERSION_MAJOR 1) -set(LASP_VERSION_MINOR 0) - -# To allow linking to static libs from other directories -cmake_policy(SET CMP0079 NEW) - -set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake") -include("BuildType") -include("QueryPythonForPybind11") - - -# Find the pybind11 package -find_pybind11_python_first() - -# This is used for code completion in vim -set(CMAKE_EXPORT_COMPILE_COMMANDS ON) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED) @@ -24,6 +8,29 @@ set(CMAKE_CXX_STANDARD_REQUIRED) option(LASP_DOUBLE_PRECISION "Compile as double precision floating point" ON) option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ON) option(LASP_HAS_ULDAQ "Compile with UlDaq backend" ON) +option(LASP_BUILD_TUNED "Tune build for current machine" OFF) + +# Use ccache if available +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}") +endif() + +# To allow linking to static libs from other directories +cmake_policy(SET CMP0079 NEW) + +# This is used for code completion in vim +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake") +include("BuildType") +include("QueryPythonForPybind11") + +# Find the pybind11 package +find_pybind11_python_first() + +find_package(OpenMP REQUIRED COMPONENTS CXX) set(LASP_FFT_BACKEND "FFTW" CACHE STRING "FFT Library backend") set_property(CACHE LASP_FFT_BACKEND PROPERTY STRINGS "FFTW" "Armadillo") @@ -33,12 +40,21 @@ set(PY_FULL_VERSION ${PROJECT_VERSION}${PY_VERSION_SUFFIX}) # Required for PYBIND11 set(CMAKE_POSITION_INDEPENDENT_CODE ON) -if(CMAKE_BUILD_TYPE STREQUAL Debug) +if(${CMAKE_BUILD_TYPE} STREQUAL "Debug") set(LASP_DEBUG True) else() set(LASP_DEBUG False) endif() +# Tune for current machine +if(LASP_BUILD_TUNED) + if(CMAKE_BUILD_TYPE STREQUAL "Debug") + 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") +endif() + # link openblas set(BLA_VENDOR OpenBLAS) find_package(BLAS REQUIRED) @@ -49,18 +65,21 @@ add_definitions(-DLASP_MAX_NFFT=33554432) # 2**25 # ####################################### End of user-adjustable variables section include(OSSpecific) - -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-type-limits") +# Enable openmp +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -Wall -Wextra -Wno-type-limits") set(CMAKE_C_FLAGS_RELEASE "-O3 -flto -mfpmath=sse -march=x86-64 -mtune=native \ -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions") # ############################# End compilation flags +include_directories(/usr/lib/python3.10/site-packages/numpy/core/include) add_subdirectory(third_party/carma) if(LASP_FFT_BACKEND STREQUAL "FFTW") - find_library(FFTW_LIBRARY NAMES fftw3 fftw) - set(FFT_LIBRARIES "${FFTW_LIBRARY}") + find_library(fftw3 REQUIRED NAMES fftw fftw3) + set(LASP_FFT_LIBS "fftw3") +message("LASP_FFT_LIBS = ${LASP_FFT_LIBS}") + elseif(LASP_FFT_BACKEND STREQUAL "Armadillo") endif() add_subdirectory(src/lasp) diff --git a/examples/Test SLM.ipynb b/examples/Test SLM.ipynb index 4891864..3c62b5d 100644 --- a/examples/Test SLM.ipynb +++ b/examples/Test SLM.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 8, + "execution_count": 1, "metadata": {}, "outputs": [ { @@ -22,26 +22,20 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "2022-08-16 21:20:04+0200 \n", - "2022-08-16 21:20:04+0200 Enter fromBiquads (lasp_slm.cpp: 82)\n", - "2022-08-16 21:20:04+0200 | Enter createBandPass (lasp_slm.cpp: 62)\n", - "2022-08-16 21:20:04+0200 | | Enter SeriesBiquad (lasp_biquadbank.cpp: 12)\n", - "2022-08-16 21:20:04+0200 | | Leave SeriesBiquad (lasp_biquadbank.cpp)\n", - "2022-08-16 21:20:04+0200 | | \n", - "2022-08-16 21:20:04+0200 | | Enter SeriesBiquad (lasp_biquadbank.cpp: 12)\n", - "2022-08-16 21:20:04+0200 | | Leave SeriesBiquad (lasp_biquadbank.cpp)\n", - "2022-08-16 21:20:04+0200 | Leave createBandPass (lasp_slm.cpp)\n", - "2022-08-16 21:20:04+0200 | \n", - "2022-08-16 21:20:04+0200 | Enter SLM (lasp_slm.cpp: 37)\n", - "2022-08-16 21:20:04+0200 | Leave SLM (lasp_slm.cpp)\n", - "2022-08-16 21:20:04+0200 Leave fromBiquads (lasp_slm.cpp)\n" + "2022-09-03 17:21:56+0200 DebugTrace-cpp 2.0.0a2 (g++ 12.2.0)\n", + "2022-09-03 17:21:56+0200 \n", + "2022-09-03 17:21:56+0200 Enter SeriesBiquad (lasp_biquadbank.cpp: 12)\n", + "2022-09-03 17:21:56+0200 Leave SeriesBiquad (lasp_biquadbank.cpp)\n", + "2022-09-03 17:21:56+0200 \n", + "2022-09-03 17:21:56+0200 Enter SeriesBiquad (lasp_biquadbank.cpp: 12)\n", + "2022-09-03 17:21:56+0200 Leave SeriesBiquad (lasp_biquadbank.cpp)\n" ] } ], @@ -64,7 +58,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -78,7 +72,7 @@ " [0., 0.]])" ] }, - "execution_count": 10, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -89,32 +83,9 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "(20100, 2)\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2022-08-16 21:20:04+0200 \n", - "2022-08-16 21:20:04+0200 Enter run (lasp_slm.cpp: 151)\n", - "2022-08-16 21:20:04+0200 | Enter filter (lasp_biquadbank.cpp: 44)\n", - "2022-08-16 21:20:04+0200 | | Enter filter (lasp_biquadbank.cpp: 44)\n", - "2022-08-16 21:20:04+0200 | | Leave filter (lasp_biquadbank.cpp)\n", - "2022-08-16 21:20:04+0200 | Leave filter (lasp_biquadbank.cpp)\n", - "2022-08-16 21:20:04+0200 | N = 0ul\n", - "2022-08-16 21:20:04+0200 | N = 0ul\n", - "2022-08-16 21:20:04+0200 Leave run (lasp_slm.cpp)\n" - ] - } - ], + "outputs": [], "source": [ "x = np.zeros(20100, dtype=float)\n", "x[:] = 1\n", @@ -221,7 +192,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.5" + "version": "3.10.6" } }, "nbformat": 4, diff --git a/examples/Test SeriesBiquad.ipynb b/examples/Test SeriesBiquad.ipynb index f0f0126..c881fa3 100644 --- a/examples/Test SeriesBiquad.ipynb +++ b/examples/Test SeriesBiquad.ipynb @@ -12,16 +12,17 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "2022-08-16 15:02:57+0200 \n", - "2022-08-16 15:02:57+0200 Enter SeriesBiquad (lasp_biquadbank.cpp: 12)\n", - "2022-08-16 15:02:57+0200 Leave SeriesBiquad (lasp_biquadbank.cpp)\n" + "2022-09-03 17:21:46+0200 DebugTrace-cpp 2.0.0a2 (g++ 12.2.0)\n", + "2022-09-03 17:21:46+0200 \n", + "2022-09-03 17:21:46+0200 Enter SeriesBiquad (lasp_biquadbank.cpp: 12)\n", + "2022-09-03 17:21:46+0200 Leave SeriesBiquad (lasp_biquadbank.cpp)\n" ] } ], @@ -32,16 +33,16 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "2022-08-16 15:02:57+0200 \n", - "2022-08-16 15:02:57+0200 Enter filter (lasp_biquadbank.cpp: 41)\n", - "2022-08-16 15:02:57+0200 Leave filter (lasp_biquadbank.cpp)\n" + "2022-09-03 17:21:47+0200 \n", + "2022-09-03 17:21:47+0200 Enter filter (lasp_biquadbank.cpp: 44)\n", + "2022-09-03 17:21:47+0200 Leave filter (lasp_biquadbank.cpp)\n" ] } ], @@ -54,7 +55,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -72,7 +73,7 @@ " [0.38742049]])" ] }, - "execution_count": 8, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -91,7 +92,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -105,7 +106,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.10.6" } }, "nbformat": 4, diff --git a/examples/Test input.ipynb b/examples/Test input.ipynb index 64894a3..2ceec77 100644 --- a/examples/Test input.ipynb +++ b/examples/Test input.ipynb @@ -1,33 +1,48 @@ { "cells": [ { - "cell_type": "code", - "execution_count": null, - "id": "b0d15138", + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "!make -j -C ~/wip/mycode/lasp" + "# Test DAQ input on a device" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# !make -j -C ~/wip/mycode/lasp" ] }, { "cell_type": "code", "execution_count": null, - "id": "1787e24c", "metadata": {}, "outputs": [], "source": [ "import lasp\n", - "ds = lasp.DeviceInfo.getDeviceInfo()" + "# Get handle to stream manager\n", + "mgr = lasp.StreamMgr.getInstance()" ] }, { "cell_type": "code", "execution_count": null, - "id": "22ae99b1", "metadata": {}, "outputs": [], "source": [ + "ds = mgr.getDeviceInfo()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Search for a device\n", "for i, d in enumerate(ds):\n", " print(f'{i}: ' + d.device_name)\n", "d = ds[1]" @@ -36,30 +51,30 @@ { "cell_type": "code", "execution_count": null, - "id": "b1321c4a", "metadata": {}, "outputs": [], "source": [ + "# Check the available block sizes\n", "d.availableFramesPerBlock" ] }, { "cell_type": "code", "execution_count": null, - "id": "47385b02", "metadata": {}, "outputs": [], "source": [ + "# Create a configuration and enable some input channels\n", "config = lasp.DaqConfiguration(d)\n", "config.inchannel_config[0].enabled = True\n", "config.inchannel_config[1].enabled = True\n", + "# Choose a different number of frames per block\n", "config.framesPerBlockIndex = 4" ] }, { "cell_type": "code", "execution_count": null, - "id": "d12f84b7", "metadata": {}, "outputs": [], "source": [ @@ -70,27 +85,16 @@ { "cell_type": "code", "execution_count": null, - "id": "902ce309", "metadata": {}, "outputs": [], "source": [ - "mgr = lasp.StreamMgr.getInstance()" + "# Start a stream with a configuration\n", + "mgr.startStream(config)" ] }, { "cell_type": "code", "execution_count": null, - "id": "b209294b", - "metadata": {}, - "outputs": [], - "source": [ - "mgr.startStream(d, config)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0830ffb5", "metadata": {}, "outputs": [], "source": [ @@ -100,7 +104,6 @@ { "cell_type": "code", "execution_count": null, - "id": "a7fddc19", "metadata": {}, "outputs": [], "source": [ @@ -110,7 +113,6 @@ { "cell_type": "code", "execution_count": null, - "id": "d11c7dae", "metadata": {}, "outputs": [], "source": [ @@ -122,7 +124,6 @@ { "cell_type": "code", "execution_count": null, - "id": "0eeb2311", "metadata": {}, "outputs": [], "source": [ @@ -134,7 +135,6 @@ { "cell_type": "code", "execution_count": null, - "id": "f893b639", "metadata": { "scrolled": true }, @@ -142,6 +142,13 @@ "source": [ "del i" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -160,7 +167,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.5" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index a59bd68..1161b85 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,4 +26,3 @@ requires = [ ] build-backend = "setuptools.build_meta" - diff --git a/src/lasp/CMakeLists.txt b/src/lasp/CMakeLists.txt index 4dc1fe8..d185739 100644 --- a/src/lasp/CMakeLists.txt +++ b/src/lasp/CMakeLists.txt @@ -31,7 +31,6 @@ pybind11_add_module(lasp_cpp MODULE lasp_cpp.cpp ) target_link_libraries(lasp_cpp PRIVATE lasp_device_lib lasp_dsp_lib - carma::carma ${OpenMP_CXX_LIBRARIES} ${LASP_FFT_LIBS}) install(TARGETS lasp_cpp DESTINATION .) diff --git a/src/lasp/device/lasp_daq.h b/src/lasp/device/lasp_daq.h index dc62b77..7e4ba35 100644 --- a/src/lasp/device/lasp_daq.h +++ b/src/lasp/device/lasp_daq.h @@ -7,6 +7,9 @@ #include /** + * \defgroup device Device interfacing + * \addtogroup device + * @{ * @brief Callback of DAQ for input data. Callback should return * false for a stop request. */ @@ -18,7 +21,8 @@ using InDaqCallback = std::function; using OutDaqCallback = std::function; /** - * @brief Base cass for all DAQ instances + * @brief Base cass for all DAQ (Data Acquisition) interfaces. A DAQ can be a + * custom device, or for example a sound card. */ class Daq : public DaqConfiguration, public DeviceInfo { @@ -90,13 +94,13 @@ public: /** * @brief Start the Daq. * - * @param cb Callback function that is called with input data as argument, - * and which expects output data as a return value. If no output data is - * required, the return value of the function should be nullptr. If no input - * data is presented, the function is called with a nullptr as argument. + * @param inCallback Function that is called with input data as argument, + * and which expects output data as a return value. If no input data is + * required, this callback should be unset. + * @param outCallback Function that is called when output data is required. + * If the stream does not require output data, this callback should be unset. */ - virtual void start(InDaqCallback inCallback, - OutDaqCallback outCallback) = 0; + virtual void start(InDaqCallback inCallback, OutDaqCallback outCallback) = 0; /** * @brief Stop the Daq device. Throws an exception if the device is not @@ -199,3 +203,4 @@ public: return (neninchannels() > 0 && nenoutchannels() > 0); } }; +/** @} */ diff --git a/src/lasp/device/lasp_daqconfig.h b/src/lasp/device/lasp_daqconfig.h index 68960e7..f11fec4 100644 --- a/src/lasp/device/lasp_daqconfig.h +++ b/src/lasp/device/lasp_daqconfig.h @@ -4,6 +4,10 @@ #include #include +/** \addtogroup device + * @{ + */ + using std::string; using std::to_string; @@ -67,17 +71,30 @@ public: }; const DataTypeDescriptor dtype_desc_fl32 = { - "32-bits floating point", 4, true, - DataTypeDescriptor::DataType::dtype_fl32}; + .name = "32-bits floating point", + .sw = 4, + .is_floating = true, + .dtype = DataTypeDescriptor::DataType::dtype_fl32}; const DataTypeDescriptor dtype_desc_fl64 = { - "64-bits floating point", 8, true, - DataTypeDescriptor::DataType::dtype_fl64}; + .name = "64-bits floating point", + .sw = 8, + .is_floating = true, + .dtype = DataTypeDescriptor::DataType::dtype_fl64}; const DataTypeDescriptor dtype_desc_int8 = { - "8-bits integer", 1, false, DataTypeDescriptor::DataType::dtype_int8}; + .name = "8-bits integer", + .sw = 1, + .is_floating = false, + .dtype = DataTypeDescriptor::DataType::dtype_int8}; const DataTypeDescriptor dtype_desc_int16 = { - "16-bits integer", 2, false, DataTypeDescriptor::DataType::dtype_int16}; + .name = "16-bits integer", + .sw = 2, + .is_floating = false, + .dtype = DataTypeDescriptor::DataType::dtype_int16}; const DataTypeDescriptor dtype_desc_int32 = { - "32-bits integer", 4, false, DataTypeDescriptor::DataType::dtype_int32}; + .name = "32-bits integer", + .sw = 4, + .is_floating = false, + .dtype = DataTypeDescriptor::DataType::dtype_int32}; const std::map dtype_map = { @@ -88,6 +105,12 @@ const std::map {DataTypeDescriptor::DataType::dtype_int32, dtype_desc_int32}, }; +/** + * @brief Class that specifies API related information. An API configuration is + * part of the DAQConfiguration and used to address a certain physical device. + * For that, a specific backend needs to be compiled in. Examples of API's are + * RtAudio and UlDaq. + */ class DaqApi { public: string apiname = "Invalid API"; @@ -127,7 +150,9 @@ const DaqApi rtaudioAsioApi("RtAudio Windows ASIO", 1, class DeviceInfo; /** - * @brief Stores channel configuration data for each channel. + * @brief Stores channel configuration data for each channel. I.e. the physical + * quantity measured, sensitivity, device-specific channel settings, a channel + * name and others. */ class DaqChannel { @@ -245,7 +270,7 @@ public: * * @param deviceinfo DeviceInfo structure */ - DaqConfiguration(const DeviceInfo &DeviceInfo); + DaqConfiguration(const DeviceInfo &deviceinfo); DaqConfiguration() {} @@ -294,3 +319,6 @@ public: */ int getLowestOutChannel() const; }; +/** + * @} + */ diff --git a/src/lasp/device/lasp_daqdata.h b/src/lasp/device/lasp_daqdata.h index f77eaf3..1700aee 100644 --- a/src/lasp/device/lasp_daqdata.h +++ b/src/lasp/device/lasp_daqdata.h @@ -5,6 +5,10 @@ #include #include +/** \addtogroup device + * @{ + */ + /** * @brief Data coming from / going to DAQ. **Non-interleaved format**, which * means data in buffer is ordered by channel: _ptr[frame+channel*nframes] @@ -126,3 +130,4 @@ public: DaqData::copyToRaw(channel, reinterpret_cast(buffer)); } }; +/** @} */ diff --git a/src/lasp/device/lasp_deviceinfo.h b/src/lasp/device/lasp_deviceinfo.h index c52fb43..de3a168 100644 --- a/src/lasp/device/lasp_deviceinfo.h +++ b/src/lasp/device/lasp_deviceinfo.h @@ -3,6 +3,9 @@ #include "lasp_types.h" #include "lasp_daqconfig.h" +/** \addtogroup device + * @{ + */ /** * @brief Structure containing device info parameters @@ -123,3 +126,4 @@ public: static std::vector getDeviceInfo(); }; +/** @} */ diff --git a/src/lasp/device/lasp_streammgr.cpp b/src/lasp/device/lasp_streammgr.cpp index 9bf67f3..9b96c0e 100644 --- a/src/lasp/device/lasp_streammgr.cpp +++ b/src/lasp/device/lasp_streammgr.cpp @@ -1,7 +1,7 @@ #define DEBUGTRACE_ENABLED +#include "lasp_streammgr.h" #include "debugtrace.hpp" #include "lasp_thread.h" -#include "lasp_streammgr.h" #include #include #include @@ -9,6 +9,7 @@ using std::cerr; using std::endl; +using rte = std::runtime_error; InDataHandler::InDataHandler(StreamMgr &mgr) : _mgr(mgr) { DEBUGTRACE_ENTER; @@ -39,12 +40,40 @@ InDataHandler::~InDataHandler() { StreamMgr &StreamMgr::getInstance() { DEBUGTRACE_ENTER; - getPool(); static StreamMgr mgr; return mgr; } -StreamMgr::StreamMgr() { DEBUGTRACE_ENTER; } +StreamMgr::StreamMgr() { + DEBUGTRACE_ENTER; +#if LASP_DEBUG == 1 + _main_thread_id = std::this_thread::get_id(); +#endif + // Trigger a scan for the available devices, in the background. + rescanDAQDevices(); +} +#if LASP_DEBUG == 1 +void StreamMgr::checkRightThread() const { + assert(std::this_thread::get_id() == _main_thread_id); +} +#endif +void StreamMgr::rescanDAQDevices(std::function callback, + bool background) { + checkRightThread(); + auto &pool = getPool(); + if (background) { + pool.push_task(&StreamMgr::rescanDAQDevices_impl, this, callback); + } else { + rescanDAQDevices_impl(callback); + } +} +void StreamMgr::rescanDAQDevices_impl(std::function callback) { + std::scoped_lock lck(_devices_mtx); + _devices = DeviceInfo::getDeviceInfo(); + if (callback) { + callback(); + } +} bool StreamMgr::inCallback(const DaqData &data) { /* DEBUGTRACE_ENTER; */ @@ -91,6 +120,7 @@ template bool fillData(DaqData &data, const vd &signal) { } void StreamMgr::setSiggen(std::shared_ptr siggen) { + checkRightThread(); std::scoped_lock lck(_siggen_mtx); _siggen = siggen; @@ -138,19 +168,40 @@ bool StreamMgr::outCallback(DaqData &data) { StreamMgr::~StreamMgr() { DEBUGTRACE_ENTER; + checkRightThread(); stopAllStreams(); } void StreamMgr::stopAllStreams() { + DEBUGTRACE_ENTER; + checkRightThread(); _inputStream.reset(); _outputStream.reset(); } -void StreamMgr::startStream(const DeviceInfo &devinfo, - const DaqConfiguration &config) { +void StreamMgr::startStream(const DaqConfiguration &config) { + DEBUGTRACE_ENTER; + checkRightThread(); bool isInput = std::count_if(config.inchannel_config.cbegin(), config.inchannel_config.cend(), [](auto &i) { return i.enabled; }); + + // Find the first device that matches with the configuration + + DeviceInfo devinfo; + { + std::scoped_lock lck(_devices_mtx); + + auto devinfo2 = + std::find_if(_devices.cbegin(), _devices.cend(), + [&config](const DeviceInfo& d) { return config.match(d); }); + if (devinfo2 == std::cend(_devices)) { + throw rte("Could not find a device with name " + config.device_name + + " in list of devices."); + } + devinfo = *devinfo2; + } + isInput |= config.monitorOutput && devinfo.hasInternalOutputMonitor; bool isOutput = std::count_if(config.outchannel_config.cbegin(), @@ -160,17 +211,16 @@ void StreamMgr::startStream(const DeviceInfo &devinfo, bool isDuplex = isInput && isOutput; if (!isInput && !isOutput) { - throw std::runtime_error("Neither input, nor output channels enabled for " - "stream. Cannotr start."); + throw rte("Neither input, nor output channels enabled for " + "stream. Cannotr start."); } if ((isDuplex || isInput) && _inputStream) { - throw std::runtime_error( - "Error: an input stream is already running. Please " - "first stop existing stream"); + throw rte("Error: an input stream is already running. Please " + "first stop existing stream"); } else if (isOutput && _outputStream) { - throw std::runtime_error("Error: output stream is already running. Please " - "first stop existing stream"); + throw rte("Error: output stream is already running. Please " + "first stop existing stream"); } InDaqCallback inCallback; @@ -202,10 +252,12 @@ void StreamMgr::startStream(const DeviceInfo &devinfo, } } void StreamMgr::stopStream(const StreamType t) { + DEBUGTRACE_ENTER; + checkRightThread(); switch (t) { case (StreamType::input): { if (!_inputStream) { - throw std::runtime_error("Input stream is not running"); + throw rte("Input stream is not running"); } _inputStream = nullptr; } break; @@ -214,28 +266,32 @@ void StreamMgr::stopStream(const StreamType t) { _inputStream = nullptr; } else { if (!_outputStream) { - throw std::runtime_error("Output stream is not running"); + throw rte("Output stream is not running"); } _outputStream = nullptr; } // end else } break; default: - throw std::runtime_error("BUG"); + throw rte("BUG"); break; } } void StreamMgr::addInDataHandler(InDataHandler &handler) { + checkRightThread(); std::scoped_lock lck(_inDataHandler_mtx); _inDataHandlers.push_back(&handler); } void StreamMgr::removeInDataHandler(InDataHandler &handler) { + checkRightThread(); std::scoped_lock lck(_inDataHandler_mtx); _inDataHandlers.remove(&handler); } Daq::StreamStatus StreamMgr::getStreamStatus(const StreamType type) const { + DEBUGTRACE_ENTER; + checkRightThread(); // Default constructor, says stream is not running, but also no errors Daq::StreamStatus s; diff --git a/src/lasp/device/lasp_streammgr.h b/src/lasp/device/lasp_streammgr.h index 5cbfa24..b182445 100644 --- a/src/lasp/device/lasp_streammgr.h +++ b/src/lasp/device/lasp_streammgr.h @@ -4,7 +4,11 @@ #include #include #include +#include +/** \addtogroup device + * @{ + */ class StreamMgr; class InDataHandler { @@ -12,7 +16,7 @@ class InDataHandler { protected: StreamMgr &_mgr; #if LASP_DEBUG == 1 - std::atomic stopCalled {false}; + std::atomic stopCalled{false}; #endif public: @@ -51,15 +55,20 @@ public: * like "pure virtual function called", or other**. */ void stop(); - }; /** * @brief Stream manager. Used to manage the input and output streams. * Implemented as a singleton: only one stream manager can be in existance for - * a given program. + * a given program. The thread that instantiates a stream manager is the only + * thread that should interact with the stream manager. If this is not true, + * undefined behavior can be expected. If LASP is compiled in debug mode, this + * fact is asserted. */ class StreamMgr { +#if LASP_DEBUG == 1 + std::thread::id _main_thread_id; +#endif /** * @brief Storage for streams. @@ -81,6 +90,9 @@ class StreamMgr { std::shared_ptr _siggen; std::mutex _siggen_mtx; + std::mutex _devices_mtx; + std::vector _devices; + StreamMgr(); friend class InDataHandler; @@ -90,7 +102,6 @@ class StreamMgr { ~StreamMgr(); public: - enum class StreamType : us { /** * @brief Input stream @@ -113,12 +124,34 @@ public: static StreamMgr &getInstance(); /** - * @brief Start a stream. + * @brief Obtain a list of devices currently available. When the StreamMgr is + * first created, this * - * @param devinfo Device information of device - * @param config Configuration of device + * @return A copy of the internal stored list of devices */ - void startStream(const DeviceInfo &devinfo, const DaqConfiguration &config); + std::vector getDeviceInfo() const { + std::scoped_lock lck(const_cast(_devices_mtx)); + return _devices; + } + + /** + * @brief Triggers a background scan of the DAQ devices, which updates the + * internally stored list of devices. + * + * @param callback Function to call when complete. + * @param background Perform searching for DAQ devices in the background. If + * set to true, the function returns immediately. + */ + void + rescanDAQDevices(std::function callback = std::function(), + bool background = false); + + /** + * @brief Start a stream based on given configuration. + * + * @param config Configuration of a device + */ + void startStream(const DaqConfiguration &config); /** * @brief Check if a certain stream is running. If running with no errors, it @@ -187,4 +220,19 @@ private: * @param handler The handler to add. */ void addInDataHandler(InDataHandler &handler); + + /** + * @brief Do the actual rescanning. + * + * @param callback + */ + void rescanDAQDevices_impl(std::function callback); + +#if LASP_DEBUG == 1 + void checkRightThread() const; +#else + void checkRightThread() const {} +#endif }; + +/** @} */ diff --git a/src/lasp/dsp/lasp_avpowerspectra.cpp b/src/lasp/dsp/lasp_avpowerspectra.cpp index 85471de..a0b6c32 100644 --- a/src/lasp/dsp/lasp_avpowerspectra.cpp +++ b/src/lasp/dsp/lasp_avpowerspectra.cpp @@ -12,10 +12,10 @@ PowerSpectra::PowerSpectra(const us nfft, const Window::WindowType w) PowerSpectra::PowerSpectra(const vd &window) : nfft(window.size()), _fft(nfft), _window(window) { - d win_pow = arma::sum(window % window); + d win_pow = arma::sum(window % window) / window.size(); /* Scale fft such that power is easily computed */ - _scale_fac = sqrt(2 / win_pow) / nfft; + _scale_fac = 2 / (win_pow * nfft * nfft); } arma::Cube PowerSpectra::compute(const dmat &input) { @@ -25,12 +25,17 @@ arma::Cube PowerSpectra::compute(const dmat &input) { // Multiply each column of the inputs element-wise with the window. input_tmp.each_col() %= _window; - cmat rfft = _fft.fft(input_tmp) * _scale_fac; + cmat rfft = _fft.fft(input_tmp); arma::cx_cube output(rfft.n_rows, input.n_cols, input.n_cols); +#pragma omp parallel for for (us i = 0; i < input.n_cols; i++) { for (us j = 0; j < input.n_cols; j++) { - output.slice(j).col(i) = rfft.col(i) % arma::conj(rfft.col(j)); + output.slice(j).col(i) = + _scale_fac * (rfft.col(i) % arma::conj(rfft.col(j))); + + output.slice(j).col(i)(0) /= 2; + output.slice(j).col(i)(nfft/2) /= 2; } } return output; diff --git a/src/lasp/dsp/lasp_avpowerspectra.h b/src/lasp/dsp/lasp_avpowerspectra.h index 558ab0f..4537e7d 100644 --- a/src/lasp/dsp/lasp_avpowerspectra.h +++ b/src/lasp/dsp/lasp_avpowerspectra.h @@ -6,6 +6,9 @@ #include #include +/** \defgroup dsp Digital Signal Processing utilities + * \ingroup dsp + * @{ /** * @brief Computes single-sided cross-power spectra for a group of channels */ @@ -151,3 +154,4 @@ public: */ us exactOverlapSamples() const { return _ps.nfft - _overlap_keep; } }; +/** @} */ diff --git a/src/lasp/dsp/lasp_biquadbank.cpp b/src/lasp/dsp/lasp_biquadbank.cpp index f3193d8..e84350e 100644 --- a/src/lasp/dsp/lasp_biquadbank.cpp +++ b/src/lasp/dsp/lasp_biquadbank.cpp @@ -6,15 +6,19 @@ using std::cerr; using std::endl; -using std::runtime_error; +using rte = std::runtime_error; SeriesBiquad::SeriesBiquad(const vd &filter_coefs) { DEBUGTRACE_ENTER; + if (filter_coefs.n_cols != 1) { + throw rte("Expected filter coefficients for a single SeriesBiquad as a " + "single column with length 6 x n_filters"); + } if (filter_coefs.n_rows % 6 != 0) { cerr << "Number of rows given: " << filter_coefs.n_rows << endl; - throw runtime_error("filter_coefs should be multiple of 6, given: " + - std::to_string(filter_coefs.n_rows)); + throw rte("filter_coefs should be multiple of 6, given: " + + std::to_string(filter_coefs.n_rows)); } us nfilters = filter_coefs.n_rows / 6; @@ -31,7 +35,7 @@ SeriesBiquad::SeriesBiquad(const vd &filter_coefs) { "absdiff", 1e-9)) { std::cerr << "Read row: " << sos.row(3) << endl; - throw std::runtime_error( + throw rte( "Filter coefficients should have fourth element (a0) equal to 1.0"); } } diff --git a/src/lasp/dsp/lasp_fft.cpp b/src/lasp/dsp/lasp_fft.cpp index 46bdef9..74e90df 100644 --- a/src/lasp/dsp/lasp_fft.cpp +++ b/src/lasp/dsp/lasp_fft.cpp @@ -7,27 +7,28 @@ // ////////////////////////////////////////////////////////////////////// #include +#include #define DEBUGTRACE_ENABLED #include "lasp_fft.h" #include "debugtrace.hpp" #include "lasp_config.h" -using std::runtime_error; +using rte = std::runtime_error; #if LASP_FFT_BACKEND == Armadillo #include "fftpack.h" class Fft_impl { -public: - us nfft; - Fft_impl(const us nfft) : nfft(nfft) { - throw runtime_error( - "This code does not output correct results, as it computes the full " - "FFT. It needs to be reworked to single-sided spectra."); + public: + us nfft; + Fft_impl(const us nfft) : nfft(nfft) { + throw runtime_error( + "This code does not output correct results, as it computes the full " + "FFT. It needs to be reworked to single-sided spectra."); - DEBUGTRACE_ENTER; - } + DEBUGTRACE_ENTER; + } - vd ifft(const vc &f) { return arma::ifft(f); } - vc fft(const vd &time) { return arma::fft(time); } + vd ifft(const vc &f) { return arma::ifft(f); } + vc fft(const vd &time) { return arma::fft(time); } }; #elif LASP_FFT_BACKEND == FFTW #include @@ -36,99 +37,120 @@ public: * @brief FFTW implementation */ class Fft_impl { -public: - us nfft; - fftw_plan forward_plan = nullptr; - fftw_plan reverse_plan = nullptr; - fftw_complex *frequencyDomain = nullptr; - d *timeDomain = nullptr; + public: + us nfft; + fftw_plan forward_plan = nullptr; + fftw_plan reverse_plan = nullptr; + fftw_complex *frequencyDomain = nullptr; + d *timeDomain = nullptr; - Fft_impl(const us nfft) : nfft(nfft) { + Fft_impl(const us nfft) : nfft(nfft) { - timeDomain = (d *)fftw_malloc(sizeof(d) * nfft); + timeDomain = (d *)fftw_malloc(sizeof(d) * nfft); - frequencyDomain = + frequencyDomain = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (nfft / 2 + 1)); - forward_plan = + forward_plan = fftw_plan_dft_r2c_1d(nfft, timeDomain, frequencyDomain, FFTW_MEASURE); - reverse_plan = + reverse_plan = fftw_plan_dft_c2r_1d(nfft, frequencyDomain, timeDomain, FFTW_MEASURE); - }; + if(!forward_plan || !reverse_plan || !timeDomain || !frequencyDomain) { + throw rte("Error allocating FFT"); + } - ~Fft_impl() { - fftw_destroy_plan(forward_plan); - fftw_destroy_plan(reverse_plan); - fftw_free(frequencyDomain); - fftw_free(timeDomain); - } + }; - vc fft(const vd &time) { + ~Fft_impl() { + fftw_destroy_plan(forward_plan); + fftw_destroy_plan(reverse_plan); + fftw_free(frequencyDomain); + fftw_free(timeDomain); + } - vc res(nfft / 2 + 1); - memcpy(timeDomain, time.memptr(), sizeof(d) * nfft); + vc fft(const vd &time) { + if(time.n_rows != nfft) { + throw rte("Invalid size of input vector, should be equal to nfft"); + } - fftw_execute(forward_plan); + vc res(nfft / 2 + 1); + memcpy(timeDomain, time.memptr(), sizeof(d) * nfft); - memcpy(reinterpret_cast(res.memptr()), frequencyDomain, - sizeof(d) * (nfft / 2 + 1)); + fftw_execute(forward_plan); - return res; - } - vd ifft(const vc &f) { - vd res(nfft); - memcpy(frequencyDomain, - reinterpret_cast(const_cast(f.memptr())), - (nfft / 2 + 1) * sizeof(c)); + memcpy(reinterpret_cast(res.memptr()), frequencyDomain, + sizeof(c) * (nfft / 2 + 1)); - fftw_execute(reverse_plan); + return res; + } + vd ifft(const vc &f) { + DEBUGTRACE_ENTER; + if(f.n_rows != nfft/2+1) { + throw rte("Invalid size of input vector, should be equal to nfft/2+1"); + } + memcpy(frequencyDomain, + reinterpret_cast(const_cast(f.memptr())), + (nfft / 2 + 1) * sizeof(c)); - memcpy(res.memptr(), frequencyDomain, nfft * sizeof(d)); + fftw_execute(reverse_plan); - /* Scale by dividing by nfft. Checked with numpy implementation - * that this indeed needs to be done for FFTW. */ - res *= 1 / nfft; + vd res(nfft); + memcpy(res.memptr(), timeDomain, nfft * sizeof(d)); - return res; - } + /* Scale by dividing by nfft. Checked with numpy implementation + * that this indeed needs to be done for FFTW. */ + res *= 1.0 / nfft; + + return res; + } }; #else #error \ - "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW" + "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW" #endif Fft::Fft(const us nfft) { DEBUGTRACE_ENTER; if (nfft == 0) { - throw std::runtime_error("Invalid nfft: 0"); + throw rte("Invalid nfft: 0"); } if (nfft >= LASP_MAX_NFFT) { - throw std::runtime_error("Invalid nfft, should be smaller than: " + - std::to_string(LASP_MAX_NFFT)); + throw rte("Invalid nfft, should be smaller than: " + + std::to_string(LASP_MAX_NFFT)); } _impl = std::make_unique(nfft); } Fft::~Fft() { } -us Fft::nfft() const { return _impl->nfft; } +us Fft::nfft() const { assert(_impl); return _impl->nfft; } + +vc Fft::fft(const vd &timedata) { assert(_impl); return _impl->fft(timedata); } + +cmat Fft::fft(const dmat &freqdata) { + DEBUGTRACE_ENTER; + assert(_impl); + cmat res(_impl->nfft/2+1, freqdata.n_cols); +#pragma omp parallel for + for (us colno = 0; colno < freqdata.n_cols; colno++) { + res.col(colno) = _impl->fft(freqdata.col(colno)); + } + return res; +} vd Fft::ifft(const vc &freqdata) { DEBUGTRACE_ENTER; - if (freqdata.size() != (_impl->nfft / 2 + 1)) { - throw runtime_error("Invalid size for frequency domain data rows." - " Should be equal to floor(nfft/2)+1"); - } + assert(_impl); return _impl->ifft(freqdata); } dmat Fft::ifft(const cmat &freqdata) { dmat res(_impl->nfft, freqdata.n_cols); +#pragma omp parallel for for (us colno = 0; colno < freqdata.n_cols; colno++) { res.col(colno) = _impl->ifft(freqdata.col(colno)); } return res; } -vc Fft::fft(const vd &timedata) { return _impl->fft(timedata); } void load_fft_wisdom(const char *wisdom) { #if LASP_FFT_BACKEND == Armadillo diff --git a/src/lasp/dsp/lasp_fft.h b/src/lasp/dsp/lasp_fft.h index 7d99489..6fe397b 100644 --- a/src/lasp/dsp/lasp_fft.h +++ b/src/lasp/dsp/lasp_fft.h @@ -28,7 +28,8 @@ char* store_fft_wisdom(void); class Fft_impl; /** - * Perform forward FFT's on real time data. Computes single-sided spectra. + * Perform forward FFT's on real time data. Computes single-sided spectra, + * equivalent to Numpy's rfft and irfft functions. */ class Fft { std::unique_ptr _impl; @@ -48,7 +49,7 @@ class Fft { /** * @brief Return nfft * - * @return nfft + * @return nfft NFFT (lenght of the DFT transform) */ us nfft() const; @@ -56,7 +57,7 @@ class Fft { * Compute the fft for a single channel of data. * * @param[in] timedata Input time data, should have size nfft - * @returns result Result complex vector + * @return Result complex vector * */ vc fft(const vd& timedata); @@ -82,8 +83,8 @@ class Fft { /** * Perform inverse FFT * - * @param[in] freqdata Frequency domain data - * @param[out] timedata Time domain result + * @param freqdata Frequency domain data + * @return timedata Time domain result */ dmat ifft(const cmat& freqdata); diff --git a/src/lasp/dsp/lasp_siggen_impl.h b/src/lasp/dsp/lasp_siggen_impl.h index 8b23478..21817e4 100644 --- a/src/lasp/dsp/lasp_siggen_impl.h +++ b/src/lasp/dsp/lasp_siggen_impl.h @@ -28,8 +28,7 @@ class Sine : public Siggen { /** * @brief Create a sine wave generator * - * @param freq_Hz - * @param level_dB + * @param freq_Hz The frequency in Hz */ Sine(const d freq_Hz); ~Sine() = default; diff --git a/src/lasp/dsp/lasp_slm.cpp b/src/lasp/dsp/lasp_slm.cpp index 3a630ed..d35014f 100644 --- a/src/lasp/dsp/lasp_slm.cpp +++ b/src/lasp/dsp/lasp_slm.cpp @@ -1,6 +1,6 @@ #define DEBUGTRACE_ENABLED -#include "lasp_slm.h" #include "debugtrace.hpp" +#include "lasp_slm.h" #include "lasp_thread.h" #include #include @@ -87,13 +87,7 @@ SLM SLM::fromBiquads(const d fs, const d Lref, const us downsampling_fac, ); } -vd SLM::run_single(const us i, d *_work, const us size) { - - vd work(_work, // Aux Memory pointer - size, // Number of elements - false, // Copy aux mem - true // Stric, means we keep being bound to the memory - ); +vd SLM::run_single(vd work,const us i) { // Filter input in-place _bandpass[i]->filter(work); @@ -156,35 +150,22 @@ dmat SLM::run(const vd &input_orig) { _pre_filter->filter(input); } - std::vector> futs; - auto &pool = getPool(); - // Fan out over multiple threads, as it is typically a heavy load dmat res(input.n_rows, _bandpass.size()); // Perform operations in-place. + +#pragma omp parallel for for (us i = 0; i < _bandpass.size(); i++) { - res.col(i) = input; - - /// It is not possible to forward a vector of values from this array in a - /// sensible way to a different thread. We therefore break it down to - /// good-old pointers and values. - futs.emplace_back( - pool.submit(&SLM::run_single, this, i, res.colptr(i), res.n_rows)); - } - + res.col(i) = run_single(input, i); /* DEBUGTRACE_PRINT(_bandpass.size()); */ /* DEBUGTRACE_PRINT(res.n_cols); */ /* DEBUGTRACE_PRINT(res.n_rows); */ /* DEBUGTRACE_PRINT(futs.size()); */ - // Wait for all threads to complete - for (us i = 0; i < _bandpass.size(); i++) { - futs[i].get(); - } - // Update the total number of samples harvested so far. NOTE: *This should be // done AFTER the threads are done!!!* + } N += input.n_rows; // Downsample, if applicable diff --git a/src/lasp/dsp/lasp_slm.h b/src/lasp/dsp/lasp_slm.h index 2a57ec3..d49fbcd 100644 --- a/src/lasp/dsp/lasp_slm.h +++ b/src/lasp/dsp/lasp_slm.h @@ -126,5 +126,5 @@ public: vd Lmax() const { return 10*arma::log10(Pmax/Lref);}; private: - vd run_single(const us i,d* inout,const us size); + vd run_single(vd input, const us filter_no); }; diff --git a/src/lasp/dsp/lasp_threadedindatahandler.h b/src/lasp/dsp/lasp_threadedindatahandler.h index 44a0603..9556ad9 100644 --- a/src/lasp/dsp/lasp_threadedindatahandler.h +++ b/src/lasp/dsp/lasp_threadedindatahandler.h @@ -4,6 +4,10 @@ const us RINGBUFFER_SIZE = 1024; +/** + * @brief Threaded in data handler. Buffers inCallback data and calls a + * callback with the same signature on a different thread. + */ class ThreadedInDataHandler: public InDataHandler { /** * @brief The queue used to push elements to the handling thread. diff --git a/src/lasp/dsp/lasp_timebuffer.cpp b/src/lasp/dsp/lasp_timebuffer.cpp index ac1bfc2..e9d515e 100644 --- a/src/lasp/dsp/lasp_timebuffer.cpp +++ b/src/lasp/dsp/lasp_timebuffer.cpp @@ -1,3 +1,5 @@ +/* #define DEBUGTRACE_ENABLED */ +#include "debugtrace.hpp" #include "lasp_timebuffer.h" #include #include @@ -16,6 +18,7 @@ class TimeBufferImp { public: void reset() { + DEBUGTRACE_ENTER; _storage.clear(); } void push(const dmat &mat) { @@ -32,6 +35,9 @@ public: } std::optional pop(const us nsamples, us keep) { + + DEBUGTRACE_ENTER; + if (keep > nsamples) throw runtime_error("keep should be <= nsamples"); diff --git a/src/lasp/dsp/lasp_types.h b/src/lasp/dsp/lasp_types.h index c350878..692a9f5 100644 --- a/src/lasp/dsp/lasp_types.h +++ b/src/lasp/dsp/lasp_types.h @@ -7,8 +7,6 @@ // needed. ////////////////////////////////////////////////////////////////////// #pragma once -#ifndef LASP_TYPES_H -#define LASP_TYPES_H #include "lasp_config.h" #include @@ -26,31 +24,30 @@ #include // true, false #include #endif + typedef size_t us; /* Size type I always use */ // To change the whole code to 32-bit floating points, change this to // float. #if LASP_FLOAT_SIZE == 32 -typedef float d; /* Shortcut for double */ +typedef float d; /* Shortcut for floating point - single precision */ #elif LASP_FLOAT_SIZE == 64 -typedef double d; /* Shortcut for double */ +typedef double d; /* Shortcut for floating point - double precision */ #else #error LASP_FLOAT_SIZE should be either 32 or 64 #endif #ifdef __cplusplus -#include -typedef std::complex c; + #include + typedef std::complex c; #else -#include -#if LASP_FLOAT_SIZE == 32 -typedef float complex c; -#else -typedef double complex c; + #include + #if LASP_FLOAT_SIZE == 32 + typedef float complex c; + #else + typedef double complex c; #endif + #endif -#endif // LASP_TYPES_H -////////////////////////////////////////////////////////////////////// - diff --git a/src/lasp/lasp_config.h.in b/src/lasp/lasp_config.h.in index c472112..f89947c 100644 --- a/src/lasp/lasp_config.h.in +++ b/src/lasp/lasp_config.h.in @@ -9,8 +9,8 @@ #ifndef LASP_CONFIG_H #define LASP_CONFIG_H -const int LASP_VERSION_MAJOR = @LASP_VERSION_MAJOR@; -const int LASP_VERSION_MINOR = @LASP_VERSION_MINOR@; +const int LASP_VERSION_MAJOR = @CMAKE_PROJECT_VERSION_MAJOR@; +const int LASP_VERSION_MINOR = @CMAKE_PROJECT_VERSION_MINOR@; /* Debug flag */ #cmakedefine01 LASP_DEBUG diff --git a/src/lasp/lasp_config.py b/src/lasp/lasp_config.py index cd5dbff..46e7af8 100644 --- a/src/lasp/lasp_config.py +++ b/src/lasp/lasp_config.py @@ -11,12 +11,12 @@ LASP_NUMPY_FLOAT_TYPE = np.float64 def zeros(shape): - return np.zeros(shape, dtype=LASP_NUMPY_FLOAT_TYPE) + return np.zeros(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order='F') def ones(shape): - return np.ones(shape, dtype=LASP_NUMPY_FLOAT_TYPE) + return np.ones(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order='F') def empty(shape): - return np.empty(shape, dtype=LASP_NUMPY_FLOAT_TYPE) + return np.empty(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order='F') diff --git a/src/lasp/lasp_slm.py b/src/lasp/lasp_slm.py index 14642c5..67b6993 100644 --- a/src/lasp/lasp_slm.py +++ b/src/lasp/lasp_slm.py @@ -4,7 +4,7 @@ Sound level meter implementation @author: J.A. de Jong - ASCEE """ -from .wrappers import Slm as pyxSlm +from .lasp_cpp import SLM as cppSLM import numpy as np from .lasp_common import (TimeWeighting, FreqWeighting, P_REF) from .filter import SPLFilterDesigner @@ -117,7 +117,7 @@ class SLM: sos = sos_overall[np.newaxis,:] self.nom_txt.append('overall') - self.slm = pyxSlm(prefilter, sos, + self.slm = cppSLM(fs, prefilter, sos, fs, tw[0], level_ref_value) dsfac = self.slm.downsampling_fac diff --git a/src/lasp/pybind11/lasp_deviceinfo.cpp b/src/lasp/pybind11/lasp_deviceinfo.cpp index 65dfc5b..a794e29 100644 --- a/src/lasp/pybind11/lasp_deviceinfo.cpp +++ b/src/lasp/pybind11/lasp_deviceinfo.cpp @@ -35,7 +35,5 @@ void init_deviceinfo(py::module& m) { devinfo.def_readonly("hasInputACCouplingSwitch", &DeviceInfo::hasInputACCouplingSwitch); - devinfo.def_static("getDeviceInfo", DeviceInfo::getDeviceInfo); - } diff --git a/src/lasp/pybind11/lasp_dsp_pybind.cpp b/src/lasp/pybind11/lasp_dsp_pybind.cpp index 3ff4002..59eae4d 100644 --- a/src/lasp/pybind11/lasp_dsp_pybind.cpp +++ b/src/lasp/pybind11/lasp_dsp_pybind.cpp @@ -1,9 +1,11 @@ -#include "carma" +#include #include "lasp_biquadbank.h" #include "lasp_siggen.h" #include "lasp_siggen_impl.h" #include "lasp_slm.h" +#include "lasp_avpowerspectra.h" #include "lasp_window.h" +#include "lasp_fft.h" #include #include @@ -12,6 +14,14 @@ namespace py = pybind11; void init_dsp(py::module &m) { + py::class_ fft(m, "Fft"); + fft.def(py::init()); + fft.def("fft", py::overload_cast(&Fft::fft)); + fft.def("fft", py::overload_cast(&Fft::fft)); + + fft.def("ifft", py::overload_cast(&Fft::ifft)); + fft.def("ifft", py::overload_cast(&Fft::ifft)); + py::class_ w(m, "Window"); py::enum_(w, "WindowType") @@ -37,7 +47,17 @@ void init_dsp(py::module &m) { return res; }); + py::class_ ps(m, "PowerSpectra"); + ps.def(py::init()); + ps.def("compute", &PowerSpectra::compute); + + py::class_ aps(m, "AvPowerSpectra"); + /* ps.def(py::init()); */ + /* ps.def("compute", &PowerSpectra::compute); */ + + py::class_ slm(m, "SLM"); + slm.def_static( "fromBiquads", py::overload_cast( diff --git a/src/lasp/pybind11/lasp_streammgr.cpp b/src/lasp/pybind11/lasp_streammgr.cpp index cc20229..bcfb7c8 100644 --- a/src/lasp/pybind11/lasp_streammgr.cpp +++ b/src/lasp/pybind11/lasp_streammgr.cpp @@ -25,4 +25,5 @@ void init_streammgr(py::module &m) { smgr.def("stopAllStreams", &StreamMgr::stopAllStreams); smgr.def("setSiggen", &StreamMgr::setSiggen); + smgr.def("getDeviceInfo", &StreamMgr::getDeviceInfo); } diff --git a/test/test_biquadbank.py b/test/test_biquadbank.py new file mode 100644 index 0000000..3809f5e --- /dev/null +++ b/test/test_biquadbank.py @@ -0,0 +1,40 @@ +#!/usr/bin/python3 +import numpy as np +from lasp import SeriesBiquad, AvPowerSpectra +from lasp.filter import SPLFilterDesigner + +import matplotlib.pyplot as plt +from scipy.signal import sosfreqz +# plt.close('all') + +# def test_cppslm2(): +# """ +# Generate a sine wave, now A-weighted +# """ +fs = 48000 +omg = 2*np.pi*1000 + +filt = SPLFilterDesigner(fs).A_Sos_design() + +bq = SeriesBiquad(filt.flatten()) + +tend=10 +t = np.linspace(0, int((tend*fs)//fs), int(tend*fs), endpoint=False) + +in_ = np.random.randn(t.size) +out = bq.filter(in_) + +from scipy.signal import welch +nfft = 48000 +freq, H1 = welch(out[:,0], + # scaling + fs=fs,nfft=nfft) + +freq, H2 = welch(in_, + # scaling + fs=fs,nfft=nfft) +# plt.figure() +plt.semilogx(freq,10*np.log10(np.abs(H1/H2))) + +omg, H_ex = sosfreqz(filt) +plt.semilogx(omg/(2*np.pi)*fs, 20*np.log10(np.abs(H_ex))) \ No newline at end of file diff --git a/test/test_cppslm.py b/test/test_cppslm.py new file mode 100644 index 0000000..f393427 --- /dev/null +++ b/test/test_cppslm.py @@ -0,0 +1,75 @@ +#!/usr/bin/python3 +import numpy as np +from lasp import SLM +from lasp.filter import SPLFilterDesigner +import matplotlib.pyplot as plt + +def test_cppslm1(): + """ + Generate a sine wave + """ + fs = 48000 + omg = 2*np.pi*1000 + + slm = SLM.fromBiquads(fs, 2e-5, 1, 0.125, [1.,0,0,1,0,0]) + + t = np.linspace(0, 10, 10*fs, endpoint=False) + + # Input signal with an rms of 1 Pa + in_ = np.sin(omg*t)*np.sqrt(2) + + # Compute overall RMS + rms = np.sqrt(np.sum(in_**2)/in_.size) + # Compute overall level + level = 20*np.log10(rms/2e-5) + + # Output of SLM + out = slm.run(in_) + + # Output of SLM should be close to theoretical + # level, at least for reasonable time constants + # (Fast, Slow etc) + assert(np.isclose(out[-1,0], level)) + + +def test_cppslm2(): + """ + Generate a sine wave, now A-weighted + """ + fs = 48000 + omg = 2*np.pi*1000 + + filt = SPLFilterDesigner(fs).A_Sos_design() + slm = SLM.fromBiquads(fs, 2e-5, 0, 0.125, filt.flatten(), [1.,0,0,1,0,0]) + + t = np.linspace(0, 10, 10*fs, endpoint=False) + + # Input signal with an rms of 1 Pa + in_ = np.sin(omg*t) *np.sqrt(2) + + # Compute overall RMS + rms = np.sqrt(np.sum(in_**2)/in_.size) + # Compute overall level + level = 20*np.log10(rms/2e-5) + + # Output of SLM + out = slm.run(in_) + + # Output of SLM should be close to theoretical + # level, at least for reasonable time constants + # (Fast, Slow etc) + assert np.isclose(out[-1,0], level, atol=1e-2) + + + +if __name__ == '__main__': + test_cppslm1() + test_cppslm2() + + # plt.plot(t,out[:,0]) + # plt.close('all') + # from scipy.signal import sosfreqz + # omg, H = sosfreqz(filt) + # freq = omg / (2*np.pi)*fs + # plt.figure() + # plt.semilogx(freq[1:],20*np.log10(np.abs(H[1:]))) diff --git a/test/test_fft.c b/test/test_fft.c deleted file mode 100644 index 2cb5311..0000000 --- a/test/test_fft.c +++ /dev/null @@ -1,33 +0,0 @@ -// test_bf.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// - -#include "lasp_fft.h" -#include "lasp_tracer.h" - - -int main() { - - setTracerLevel(0); - - iVARTRACE(15,getTracerLevel()); - - Fft* fft = Fft_create(100000); - - /* Fft_fft(fft,NULL,NULL); */ - - Fft_free(fft); - - return 0; -} - - - - -////////////////////////////////////////////////////////////////////// - - diff --git a/test/test_fft.py b/test/test_fft.py index 9cb248b..2f3b087 100644 --- a/test/test_fft.py +++ b/test/test_fft.py @@ -43,3 +43,6 @@ def test_backward_fft(): assert(np.isclose(np.linalg.norm(sig_py- sig_lasp), 0)) +if __name__ == '__main__': + test_forward_fft() + # test_backward_fft()() \ No newline at end of file diff --git a/test/test_math.c b/test/test_math.c deleted file mode 100644 index 404eeac..0000000 --- a/test/test_math.c +++ /dev/null @@ -1,50 +0,0 @@ -// test_bf.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// - -#include "lasp_fft.h" -#include "lasp_mat.h" -#include "lasp_tracer.h" -#include "lasp_alg.h" - -int main() { - - iVARTRACE(15,getTracerLevel()); - - c a[5]; - c_set(a,1-I,5); - - /* print_vc(&a); */ - - vc b = vc_alloc(5); - vc_set(&b,2); - - - printf("b:\n"); - print_vc(&b); - - - vc c1 = vc_alloc(5); - /* vc_set(&c1,10); */ - /* c_add_to(c1.ptr,a.ptr,1,3); */ - c_hadamard(c1._data,a,b._data,5); - - - - printf("c1:\n"); - print_vc(&c1); - - vc_free(&b); - - - - return 0; -} - -////////////////////////////////////////////////////////////////////// - - diff --git a/test/test_ps.py b/test/test_ps.py index ed7859f..fbeb57e 100644 --- a/test/test_ps.py +++ b/test/test_ps.py @@ -6,51 +6,140 @@ Created on Mon Jan 15 19:45:33 2018 @author: anne """ import numpy as np -from beamforming import Fft, PowerSpectra, cls -cls -nfft=2048 -print('nfft:',nfft) -#print(nfft) -nchannels = 2 +from lasp import PowerSpectra, Window +import matplotlib.pyplot as plt +plt.close('all') +# def test_ps(): -t = np.linspace(0,1,nfft+1) -# print(t) -x1 = (np.cos(4*np.pi*t[:-1])+3.2*np.sin(6*np.pi*t[:-1]))[:,np.newaxis]+10 -x = np.vstack([x1.T]*nchannels).T -# Using transpose to get the strides right -x = np.random.randn(nchannels,nfft).T -print("strides: ",x.strides) -# x.strides = (8,nfft*8)x -# print("signal:",x) +nfft = 8 +t = np.linspace(0, 1.0, nfft, endpoint=False) -xms = np.sum(x**2,axis=0)/nfft -print('Total signal power time domain: ', xms) +ps = PowerSpectra(nfft, Window.Rectangular) -X = np.fft.rfft(x,axis=0) -# X =np.fft.fft(x) -#X =np.fft.rfft(x) - -# print(X) -Xs = 2*X/nfft -Xs[np.where(np.abs(Xs) < 1e-10)] = 0 -Xs[0] /= np.sqrt(2) -Xs[-1] /= np.sqrt(2) -# print('single sided amplitude spectrum:\n',Xs) +sig = np.random.randn(nfft) +freq = 4 +omg = 2*np.pi*freq +# sig = 8*np.cos(omg*t) -power = Xs*np.conj(Xs)/2 +cps = ps.compute(sig) -# print('Frequency domain signal power\n', power) -print('Total signal power', np.sum(power,axis=0).real) +pow1 = np.sum(sig**2)/sig.size +pow2 = np.sum((cps[:,0,0]).real) -pstest = PowerSpectra(nfft,nchannels) -ps = pstest.compute(x) +# print(pow1) +# print(pow2) +plt.plot(cps[:,0,0]) +assert np.isclose(pow2 - pow1,0, atol=1e-1) +# test_ps() +# plt.plot(res_lasp.real-res_npy.real) +# plt.plot(res_lasp.imag-res_npy.imag) +# plt.plot(res_npy.real) +# plt.plot(res_npy.imag) + # plt.plot(t, sig) +# print('nfft:',nfft) +# #print(nfft) +# nchannels = 2 -fft = Fft(nfft,nchannels) -fft.fft(x) +# t = np.linspace(0,1,nfft+1) +# # print(t) +# x1 = (np.cos(4*np.pi*t[:-1])+3.2*np.sin(6*np.pi*t[:-1]))[:,np.newaxis]+10 +# x = np.vstack([x1.T]*nchannels).T +# # Using transpose to get the strides right +# x = np.random.randn(nchannels,nfft).T +# print("strides: ",x.strides) +# # x.strides = (8,nfft*8)x +# # print("signal:",x) -ps[np.where(np.abs(ps) < 1e-10)] = 0+0j +# xms = np.sum(x**2,axis=0)/nfft +# print('Total signal power time domain: ', xms) -print('our ps: \n' , ps) +# X = np.fft.rfft(x,axis=0) +# # X =np.fft.fft(x) +# #X =np.fft.rfft(x) -print('Our total signal power: ',np.sum(ps,axis=0).real) +# # print(X) +# Xs = 2*X/nfft +# Xs[np.where(np.abs(Xs) < 1e-10)] = 0 +# Xs[0] /= np.sqrt(2) +# Xs[-1] /= np.sqrt(2) +# # print('single sided amplitude spectrum:\n',Xs) + + +# power = Xs*np.conj(Xs)/2 + +# # print('Frequency domain signal power\n', power) +# print('Total signal power', np.sum(power,axis=0).real) + +# pstest = PowerSpectra(nfft,nchannels) +# ps = pstest.compute(x) + +# fft = Fft(nfft,nchannels) +# fft.fft(x) + +# ps[np.where(np.abs(ps) < 1e-10)] = 0+0j + +# print('our ps: \n' , ps) + +# print('Our total signal power: ',np.sum(ps,axis=0).real) + + +# if __name__ == '__main__': +# nfft=2048 +# fs = 2048 +# ps = PowerSpectra(nfft, Window.Rectangular) + +# t = np.linspace(0, 1.0, nfft, endpoint=False) +# freq = 10 +# omg = 2*np.pi*freq +# sig = np.sin(omg*t) + +# res = ps.compute(sig) + +# plt.plot(res[:,0,0]) +# # plt.plot(t, sig) +# print('nfft:',nfft) +# #print(nfft) +# nchannels = 2 + +# t = np.linspace(0,1,nfft+1) +# # print(t) +# x1 = (np.cos(4*np.pi*t[:-1])+3.2*np.sin(6*np.pi*t[:-1]))[:,np.newaxis]+10 +# x = np.vstack([x1.T]*nchannels).T +# # Using transpose to get the strides right +# x = np.random.randn(nchannels,nfft).T +# print("strides: ",x.strides) +# # x.strides = (8,nfft*8)x +# # print("signal:",x) + +# xms = np.sum(x**2,axis=0)/nfft +# print('Total signal power time domain: ', xms) + +# X = np.fft.rfft(x,axis=0) +# # X =np.fft.fft(x) +# #X =np.fft.rfft(x) + +# # print(X) +# Xs = 2*X/nfft +# Xs[np.where(np.abs(Xs) < 1e-10)] = 0 +# Xs[0] /= np.sqrt(2) +# Xs[-1] /= np.sqrt(2) +# # print('single sided amplitude spectrum:\n',Xs) + + +# power = Xs*np.conj(Xs)/2 + +# # print('Frequency domain signal power\n', power) +# print('Total signal power', np.sum(power,axis=0).real) + +# pstest = PowerSpectra(nfft,nchannels) +# ps = pstest.compute(x) + +# fft = Fft(nfft,nchannels) +# fft.fft(x) + +# ps[np.where(np.abs(ps) < 1e-10)] = 0+0j + +# print('our ps: \n' , ps) + +# print('Our total signal power: ',np.sum(ps,axis=0).real) diff --git a/test/test_rtaudio.py b/test/test_rtaudio.py deleted file mode 100644 index 582986c..0000000 --- a/test/test_rtaudio.py +++ /dev/null @@ -1,35 +0,0 @@ -#!/usr/bin/python3 -import numpy as np -from lasp_rtaudio import RtAudio, SampleFormat, Format_SINT32, Format_FLOAT64 -import time - - -nframes = 0 -samplerate = 48000 -omg = 2*np.pi*1000 - -def mycallback(input_, nframes, streamtime): - t = np.linspace(streamtime, streamtime + nframes/samplerate, - nframes)[np.newaxis,:] - outp = 0.1*np.sin(omg*t) - return outp, 0 - -if __name__ == '__main__': - pa = RtAudio() - count = pa.getDeviceCount() - # dev = pa.getDeviceInfo(0) - for i in range(count): - dev = pa.getDeviceInfo(i) - print(dev) - - outputparams = {'deviceid': 0, 'nchannels': 1, 'firstchannel': 0} - pa.openStream(outputparams, None , Format_FLOAT64,samplerate, 512, mycallback) - pa.startStream() - - input() - - pa.stopStream() - pa.closeStream() - - -