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.
This commit is contained in:
parent
10749137ec
commit
70891ceaf4
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -26,4 +26,3 @@ requires = [
|
|||
]
|
||||
build-backend = "setuptools.build_meta"
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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 .)
|
||||
|
|
|
|||
|
|
@ -7,6 +7,9 @@
|
|||
#include <memory>
|
||||
|
||||
/**
|
||||
* \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<bool(const DaqData &)>;
|
|||
using OutDaqCallback = std::function<bool(DaqData &)>;
|
||||
|
||||
/**
|
||||
* @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);
|
||||
}
|
||||
};
|
||||
/** @} */
|
||||
|
|
|
|||
|
|
@ -4,6 +4,10 @@
|
|||
#include <map>
|
||||
#include <vector>
|
||||
|
||||
/** \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<DataTypeDescriptor::DataType, const DataTypeDescriptor>
|
||||
dtype_map = {
|
||||
|
|
@ -88,6 +105,12 @@ const std::map<DataTypeDescriptor::DataType, const DataTypeDescriptor>
|
|||
{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;
|
||||
};
|
||||
/**
|
||||
* @}
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -5,6 +5,10 @@
|
|||
#include <gsl/gsl-lite.hpp>
|
||||
#include <memory>
|
||||
|
||||
/** \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<uint8_t *>(buffer));
|
||||
}
|
||||
};
|
||||
/** @} */
|
||||
|
|
|
|||
|
|
@ -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<DeviceInfo> getDeviceInfo();
|
||||
};
|
||||
|
||||
/** @} */
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
#define DEBUGTRACE_ENABLED
|
||||
#include "lasp_streammgr.h"
|
||||
#include "debugtrace.hpp"
|
||||
#include "lasp_thread.h"
|
||||
#include "lasp_streammgr.h"
|
||||
#include <algorithm>
|
||||
#include <assert.h>
|
||||
#include <functional>
|
||||
|
|
@ -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<void()> 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<void()> 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 <typename T> bool fillData(DaqData &data, const vd &signal) {
|
|||
}
|
||||
|
||||
void StreamMgr::setSiggen(std::shared_ptr<Siggen> siggen) {
|
||||
checkRightThread();
|
||||
|
||||
std::scoped_lock<std::mutex> 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<std::mutex> lck(_inDataHandler_mtx);
|
||||
_inDataHandlers.push_back(&handler);
|
||||
}
|
||||
void StreamMgr::removeInDataHandler(InDataHandler &handler) {
|
||||
checkRightThread();
|
||||
std::scoped_lock<std::mutex> 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;
|
||||
|
||||
|
|
|
|||
|
|
@ -4,7 +4,11 @@
|
|||
#include <list>
|
||||
#include <memory>
|
||||
#include <mutex>
|
||||
#include <thread>
|
||||
|
||||
/** \addtogroup device
|
||||
* @{
|
||||
*/
|
||||
class StreamMgr;
|
||||
|
||||
class InDataHandler {
|
||||
|
|
@ -12,7 +16,7 @@ class InDataHandler {
|
|||
protected:
|
||||
StreamMgr &_mgr;
|
||||
#if LASP_DEBUG == 1
|
||||
std::atomic<bool> stopCalled {false};
|
||||
std::atomic<bool> 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> _siggen;
|
||||
std::mutex _siggen_mtx;
|
||||
|
||||
std::mutex _devices_mtx;
|
||||
std::vector<DeviceInfo> _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<DeviceInfo> getDeviceInfo() const {
|
||||
std::scoped_lock lck(const_cast<std::mutex&>(_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<void()> callback = std::function<void()>(),
|
||||
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<void()> callback);
|
||||
|
||||
#if LASP_DEBUG == 1
|
||||
void checkRightThread() const;
|
||||
#else
|
||||
void checkRightThread() const {}
|
||||
#endif
|
||||
};
|
||||
|
||||
/** @} */
|
||||
|
|
|
|||
|
|
@ -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<c> PowerSpectra::compute(const dmat &input) {
|
||||
|
|
@ -25,12 +25,17 @@ arma::Cube<c> 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;
|
||||
|
|
|
|||
|
|
@ -6,6 +6,9 @@
|
|||
#include <memory>
|
||||
#include <optional>
|
||||
|
||||
/** \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; }
|
||||
};
|
||||
/** @} */
|
||||
|
|
|
|||
|
|
@ -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");
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -7,27 +7,28 @@
|
|||
//
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
#include <memory>
|
||||
#include <cassert>
|
||||
#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 <fftw3.h>
|
||||
|
|
@ -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<void *>(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<void *>(const_cast<c *>(f.memptr())),
|
||||
(nfft / 2 + 1) * sizeof(c));
|
||||
memcpy(reinterpret_cast<void *>(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<void *>(const_cast<c *>(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<Fft_impl>(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
|
||||
|
|
|
|||
|
|
@ -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<Fft_impl> _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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
#define DEBUGTRACE_ENABLED
|
||||
#include "lasp_slm.h"
|
||||
#include "debugtrace.hpp"
|
||||
#include "lasp_slm.h"
|
||||
#include "lasp_thread.h"
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
|
|
@ -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<std::future<vd>> 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
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
};
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -1,3 +1,5 @@
|
|||
/* #define DEBUGTRACE_ENABLED */
|
||||
#include "debugtrace.hpp"
|
||||
#include "lasp_timebuffer.h"
|
||||
#include <algorithm>
|
||||
#include <cassert>
|
||||
|
|
@ -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<dmat> pop(const us nsamples, us keep) {
|
||||
|
||||
DEBUGTRACE_ENTER;
|
||||
|
||||
if (keep > nsamples)
|
||||
throw runtime_error("keep should be <= nsamples");
|
||||
|
||||
|
|
|
|||
|
|
@ -7,8 +7,6 @@
|
|||
// needed.
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
#pragma once
|
||||
#ifndef LASP_TYPES_H
|
||||
#define LASP_TYPES_H
|
||||
#include "lasp_config.h"
|
||||
#include <stddef.h>
|
||||
|
||||
|
|
@ -26,31 +24,30 @@
|
|||
#include <stdbool.h> // true, false
|
||||
#include <stddef.h>
|
||||
#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 <complex>
|
||||
typedef std::complex<d> c;
|
||||
#include <complex>
|
||||
typedef std::complex<d> c;
|
||||
#else
|
||||
#include <complex.h>
|
||||
#if LASP_FLOAT_SIZE == 32
|
||||
typedef float complex c;
|
||||
#else
|
||||
typedef double complex c;
|
||||
#include <complex.h>
|
||||
#if LASP_FLOAT_SIZE == 32
|
||||
typedef float complex c;
|
||||
#else
|
||||
typedef double complex c;
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
#endif // LASP_TYPES_H
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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')
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -35,7 +35,5 @@ void init_deviceinfo(py::module& m) {
|
|||
devinfo.def_readonly("hasInputACCouplingSwitch",
|
||||
&DeviceInfo::hasInputACCouplingSwitch);
|
||||
|
||||
devinfo.def_static("getDeviceInfo", DeviceInfo::getDeviceInfo);
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1,9 +1,11 @@
|
|||
#include "carma"
|
||||
#include <carma>
|
||||
#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 <iostream>
|
||||
#include <pybind11/pybind11.h>
|
||||
|
||||
|
|
@ -12,6 +14,14 @@ namespace py = pybind11;
|
|||
|
||||
void init_dsp(py::module &m) {
|
||||
|
||||
py::class_<Fft> fft(m, "Fft");
|
||||
fft.def(py::init<us>());
|
||||
fft.def("fft", py::overload_cast<const vd&>(&Fft::fft));
|
||||
fft.def("fft", py::overload_cast<const dmat&>(&Fft::fft));
|
||||
|
||||
fft.def("ifft", py::overload_cast<const vc&>(&Fft::ifft));
|
||||
fft.def("ifft", py::overload_cast<const cmat&>(&Fft::ifft));
|
||||
|
||||
py::class_<Window> w(m, "Window");
|
||||
|
||||
py::enum_<Window::WindowType>(w, "WindowType")
|
||||
|
|
@ -37,7 +47,17 @@ void init_dsp(py::module &m) {
|
|||
return res;
|
||||
});
|
||||
|
||||
py::class_<PowerSpectra> ps(m, "PowerSpectra");
|
||||
ps.def(py::init<us, const Window::WindowType>());
|
||||
ps.def("compute", &PowerSpectra::compute);
|
||||
|
||||
py::class_<AvPowerSpectra> aps(m, "AvPowerSpectra");
|
||||
/* ps.def(py::init<us, const Window::WindowType>()); */
|
||||
/* ps.def("compute", &PowerSpectra::compute); */
|
||||
|
||||
|
||||
py::class_<SLM> slm(m, "SLM");
|
||||
|
||||
slm.def_static(
|
||||
"fromBiquads",
|
||||
py::overload_cast<const d, const d, const us, const d, const dmat &>(
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
|
|
|||
40
test/test_biquadbank.py
Normal file
40
test/test_biquadbank.py
Normal file
|
|
@ -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)))
|
||||
75
test/test_cppslm.py
Normal file
75
test/test_cppslm.py
Normal file
|
|
@ -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:])))
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
|
|
@ -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()()
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
163
test/test_ps.py
163
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)
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user