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() - - -