From 7095f9d5e7b72d3fbf9dba98701aeb7d485a96ea Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Wed, 29 Jun 2022 12:25:32 +0200 Subject: [PATCH] Intermediate commit. Ready for some serious testing. --- .gitignore | 1 + .gitmodules | 3 + CMakeLists.txt | 65 +-- armadillo-code | 1 + lasp/CMakeLists.txt | 13 +- lasp/__init__.py | 1 - lasp/c/CMakeLists.txt | 45 +- lasp/c/lasp_fft.h | 6 + lasp/c/lasp_math_raw.h | 41 +- lasp/c/lasp_nprocs.h | 19 +- lasp/c/lasp_siggen.c | 472 ------------------ lasp/c/lasp_siggen.h | 100 ---- lasp/c/lasp_tracer.c | 78 --- lasp/c/lasp_tracer.h | 238 --------- lasp/c/lasp_window.c | 111 ----- lasp/c/lasp_window.h | 37 -- lasp/config.pxi.in | 10 - lasp/device/CMakeLists.txt | 30 +- lasp/device/__init__.py | 5 +- lasp/device/lasp_daq.cpp | 3 + lasp/device/lasp_daq.h | 38 +- lasp/device/lasp_daqconfig.h | 6 +- lasp/device/lasp_daqdata.h | 4 +- lasp/device/lasp_devicepybind.cpp | 161 +----- lasp/device/lasp_devicepybind.h | 44 -- lasp/device/lasp_rtaudiodaq.cpp | 12 +- lasp/device/lasp_streammgr.cpp | 151 ++++++ lasp/device/lasp_streammgr.h | 187 +++++++ lasp/device/lasp_uldaq.cpp | 1 - lasp/dsp/CMakeLists.txt | 21 + lasp/dsp/__init__.py | 1 + lasp/{c => dsp}/lasp_config.h.in | 0 lasp/dsp/lasp_dsp_pybind.cpp | 20 + lasp/dsp/lasp_filter.cpp | 0 lasp/dsp/lasp_filter.h | 14 + lasp/dsp/lasp_mathtypes.h | 45 ++ lasp/dsp/lasp_siggen.cpp | 42 ++ lasp/dsp/lasp_siggen.h | 77 +++ lasp/dsp/lasp_siggen_impl.cpp | 245 +++++++++ lasp/dsp/lasp_siggen_impl.h | 79 +++ lasp/{c => dsp}/lasp_types.h | 8 - lasp/dsp/lasp_window.cpp | 57 +++ lasp/dsp/lasp_window.h | 73 +++ lasp/lasp_common.py | 12 +- lasp/wrappers.pyx | 798 ------------------------------ setup.py | 0 46 files changed, 1182 insertions(+), 2193 deletions(-) create mode 160000 armadillo-code delete mode 100644 lasp/c/lasp_siggen.c delete mode 100644 lasp/c/lasp_siggen.h delete mode 100644 lasp/c/lasp_tracer.c delete mode 100644 lasp/c/lasp_tracer.h delete mode 100644 lasp/c/lasp_window.c delete mode 100644 lasp/c/lasp_window.h delete mode 100644 lasp/config.pxi.in delete mode 100644 lasp/device/lasp_devicepybind.h create mode 100644 lasp/device/lasp_streammgr.cpp create mode 100644 lasp/device/lasp_streammgr.h create mode 100644 lasp/dsp/CMakeLists.txt create mode 100644 lasp/dsp/__init__.py rename lasp/{c => dsp}/lasp_config.h.in (100%) create mode 100644 lasp/dsp/lasp_dsp_pybind.cpp create mode 100644 lasp/dsp/lasp_filter.cpp create mode 100644 lasp/dsp/lasp_filter.h create mode 100644 lasp/dsp/lasp_mathtypes.h create mode 100644 lasp/dsp/lasp_siggen.cpp create mode 100644 lasp/dsp/lasp_siggen.h create mode 100644 lasp/dsp/lasp_siggen_impl.cpp create mode 100644 lasp/dsp/lasp_siggen_impl.h rename lasp/{c => dsp}/lasp_types.h (87%) create mode 100644 lasp/dsp/lasp_window.cpp create mode 100644 lasp/dsp/lasp_window.h delete mode 100644 lasp/wrappers.pyx mode change 100644 => 100755 setup.py diff --git a/.gitignore b/.gitignore index 7fd86ff..798fb0a 100644 --- a/.gitignore +++ b/.gitignore @@ -32,3 +32,4 @@ lasp/device/lasp_daq.cxx lasp/c/lasp_config.h compile_commands.json .cache +lasp_config.h diff --git a/.gitmodules b/.gitmodules index e9c585f..14931f2 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,3 +7,6 @@ [submodule "DebugTrace-cpp"] path = DebugTrace-cpp url = https://github.com/MasatoKokubo/DebugTrace-cpp +[submodule "armadillo-code"] + path = armadillo-code + url = https://gitlab.com/conradsnicta/armadillo-code diff --git a/CMakeLists.txt b/CMakeLists.txt index f432355..ee90237 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required (VERSION 3.12) +cmake_minimum_required (VERSION 3.16) # To allow linking to static libs from other directories cmake_policy(SET CMP0079 NEW) @@ -6,28 +6,25 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/lasp/cma include("BuildType") # This is used for code completion in vim -set(CMAKE_EXPORT_COMPILE_COMMANDS=ON) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) project(LASP LANGUAGES C CXX) set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED) set(CMAKE_C_STANDARD 11) -# Whether we want to use blas yes or no -option(LASP_USE_BLAS "Use external blas library for math" ON) -option(LASP_ARM "Compile subset of code for ARM real time (Bela board)" OFF) option(LASP_DOUBLE_PRECISION "Compile as double precision floating point" ON) -option(LASP_PARALLEL "Parallel processing" ON) option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ON) option(LASP_HAS_ULDAQ "Compile with UlDaq backend" ON) -set(LASP_MAX_NUM_CHANNELS 16 CACHE STRING "Maximum number of audio channels that is compiled in in code") -set(LASP_MAX_NUM_THREADS 30 CACHE STRING "The maximum number of simultaneous threads in parallel processing") set(LASP_TRACERNAME "defaulttracer" CACHE STRING "Name of tracer variable containing level") set(LASP_FFT_BACKEND "FFTW" CACHE STRING "FFT Library backend") set_property(CACHE LASP_FFT_BACKEND PROPERTY STRINGS "FFTW" "FFTPack" "None") -set(PYBIND11_FINDPYTHON ON) +find_package(Python3 COMPONENTS Interpreter Development REQUIRED) +include_directories(${Python3_INCLUDE_DIRS}) find_package(pybind11 CONFIG REQUIRED) + # Required for PYBIND11 set(POSITION_INDEPENDENT_CODE True) @@ -44,11 +41,12 @@ else() set(LASP_THREADING_LIBRARIES "") endif() -if(LASP_USE_BLAS) - # link openblas - set(BLA_VENDOR OpenBLAS) - find_package(BLAS REQUIRED) -endif() +# link openblas +set(BLA_VENDOR OpenBLAS) +find_package(BLAS REQUIRED) + +# Armadillo +include_directories(SYSTEM armadillo-code/include) # Reasonable maximum to the nfft size, at 48kHz this is 700s of data... @@ -67,6 +65,7 @@ endif() # General make flags set(CMAKE_POSITION_INDEPENDENT_CODE ON) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wextra -Wno-type-limits \ -Werror=implicit-function-declaration -Wno-unused-parameter \ -Werror=return-type -Wfatal-errors") @@ -98,7 +97,6 @@ else() # Linux compile include_directories(/usr/local/include/rtaudio) include_directories(/usr/include/rtaudio) link_directories(/usr/local/lib) - # This should become optional later on, and be added to the windows list as # well. @@ -108,21 +106,12 @@ endif() set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-type-limits") -# Debug make flags -set(CMAKE_C_FLAGS_DEBUG "-g" ) - set(CMAKE_C_FLAGS_RELEASE "-O2 -mfpmath=sse -march=x86-64 -mtune=native \ -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions") # ############################# End compilation flags -# Python searching. -# set(Python_ADDITIONAL_VERSIONS "3.8") -# set(python_version_windll "38") -find_package(PythonLibs REQUIRED ) -find_package(PythonInterp REQUIRED) - # Add FFTpack dir if used as FFT backend if(LASP_FFTPACK_BACKEND) add_subdirectory(fftpack) @@ -131,10 +120,10 @@ if(LASP_FFTPACK_BACKEND) ) endif() -include_directories(lasp/c) include_directories(STL-Threadsafe/include) include_directories(gsl-lite/include) include_directories(DebugTrace-cpp/include) + add_subdirectory(lasp) add_subdirectory(gsl-lite) add_subdirectory(test) @@ -148,21 +137,21 @@ set(DEPS "${CMAKE_CURRENT_SOURCE_DIR}/*.py" # "lasp_device_lib") ) -set(OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp") +# set(OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp") -# configure_file(${SETUP_PY_IN} ${SETUP_PY}) -add_custom_command(OUTPUT ${OUTPUT} - COMMAND ${PYTHON_EXECUTABLE} ${SETUP_PY} build - COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT} - DEPENDS ${DEPS}) +# # configure_file(${SETUP_PY_IN} ${SETUP_PY}) +# add_custom_command(OUTPUT ${OUTPUT} +# COMMAND ${PYTHON_EXECUTABLE} ${SETUP_PY} build +# COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT} +# DEPENDS ${DEPS}) -add_custom_target(target ALL DEPENDS ${OUTPUT}) +# add_custom_target(target ALL DEPENDS ${OUTPUT}) -if(DEFINED INSTALL_DEBUG) - set(EXTRA_SETUP_ARG --user -e) -else() - set(EXTRA_SETUP_ARG "") -endif() +# if(DEFINED INSTALL_DEBUG) +# set(EXTRA_SETUP_ARG --user -e) +# else() +# set(EXTRA_SETUP_ARG "") +# endif() -install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pip install ${EXTRA_SETUP_ARG} .)") +# install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pip install ${EXTRA_SETUP_ARG} .)") diff --git a/armadillo-code b/armadillo-code new file mode 160000 index 0000000..b572743 --- /dev/null +++ b/armadillo-code @@ -0,0 +1 @@ +Subproject commit b5727433d342afca53aca0ee16ecf1caaa661821 diff --git a/lasp/CMakeLists.txt b/lasp/CMakeLists.txt index 447105e..db5719f 100644 --- a/lasp/CMakeLists.txt +++ b/lasp/CMakeLists.txt @@ -1,14 +1,3 @@ -configure_file(config.pxi.in config.pxi) - -# This is used for code completion in vim -set(CMAKE_EXPORT_COMPILE_COMMANDS ON) - add_subdirectory(c) add_subdirectory(device) - -# set_source_files_properties(wrappers.c PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} ${CYTHON_EXTRA_C_FLAGS}") -# cython_add_module(wrappers wrappers.pyx) -# target_link_libraries(wrappers lasp_lib ${LASP_THREADING_LIBRARIES}) -# if(win32) -# target_link_libraries(wrappers python${python_version_windll}) -# endif(win32) +add_subdirectory(dsp) diff --git a/lasp/__init__.py b/lasp/__init__.py index a5254ee..b879372 100644 --- a/lasp/__init__.py +++ b/lasp/__init__.py @@ -2,7 +2,6 @@ from .lasp_common import * # P_REF, FreqWeighting, TimeWeighting, getTime, getFreq, Qty, SIQtys, Window, lasp_shelve, this_lasp_shelve, W_REF, U_REF, I_REF, dBFS_REF, AvType from .lasp_avstream import * # StreamManager, ignoreSigInt, StreamStatus -from .wrappers import * # AvPowerSpectra, SosFilterBank, FilterBank, Siggen, sweep_flag_forward, sweep_flag_backward, sweep_flag_linear, sweep_flag_exponential, load_fft_wisdom, store_fft_wisdom from .lasp_atomic import * # Atomic from .lasp_imptube import * # TwoMicImpedanceTube from .lasp_measurement import * # Measurement, scaleBlockSens diff --git a/lasp/c/CMakeLists.txt b/lasp/c/CMakeLists.txt index ff86eb5..905866d 100644 --- a/lasp/c/CMakeLists.txt +++ b/lasp/c/CMakeLists.txt @@ -1,30 +1,29 @@ -configure_file(lasp_config.h.in lasp_config.h) if(!LASP_DEBUG) SET_SOURCE_FILES_PROPERTIES(lasp_sosfilterbank.c PROPERTIES COMPILE_FLAGS -O3) SET_SOURCE_FILES_PROPERTIES(lasp_slm.c PROPERTIES COMPILE_FLAGS -O3) SET_SOURCE_FILES_PROPERTIES(lasp_eq.c PROPERTIES COMPILE_FLAGS -O3) endif(!LASP_DEBUG) -add_library(lasp_lib - lasp_fft.c - lasp_mat.c - lasp_math_raw.c - lasp_alg.c - lasp_assert.c - lasp_tracer.c - lasp_window.c - lasp_aps.c - lasp_ps.c - lasp_mq.c - lasp_siggen.c - lasp_worker.c - lasp_nprocs.c - lasp_dfifo.c - lasp_firfilterbank.c - lasp_sosfilterbank.c - lasp_decimation.c - lasp_slm.c - lasp_eq.c - ) -target_link_libraries(lasp_lib ${FFT_LIBRARIES} ${BLAS_LIBRARIES}) +# add_library(lasp_lib +# lasp_fft.c +# lasp_mat.c +# # lasp_math_raw.c +# # lasp_alg.c +# lasp_assert.c +# lasp_tracer.c +# # lasp_window.c +# lasp_aps.c +# lasp_ps.c +# # lasp_mq.c +# lasp_worker.c +# lasp_nprocs.c +# lasp_dfifo.c +# lasp_firfilterbank.c +# lasp_sosfilterbank.c +# lasp_decimation.c +# lasp_slm.c +# lasp_eq.c +# ) + +# target_link_libraries(lasp_lib ${FFT_LIBRARIES} ${BLAS_LIBRARIES}) diff --git a/lasp/c/lasp_fft.h b/lasp/c/lasp_fft.h index 701f658..5b1f39d 100644 --- a/lasp/c/lasp_fft.h +++ b/lasp/c/lasp_fft.h @@ -10,6 +10,9 @@ #define LASP_FFT_H #include "lasp_types.h" #include "lasp_mat.h" +#ifdef __cplusplus +extern "C" { +#endif /** * Load wisdom data from a NULL-terminated character string. Note that at this @@ -99,5 +102,8 @@ void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata); */ void Fft_free(Fft* fft); +#ifdef __cplusplus +} +#endif #endif // LASP_FFT_H ////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_math_raw.h b/lasp/c/lasp_math_raw.h index ed91b87..d888f7b 100644 --- a/lasp/c/lasp_math_raw.h +++ b/lasp/c/lasp_math_raw.h @@ -8,6 +8,10 @@ #pragma once #ifndef LASP_MATH_RAW_H #define LASP_MATH_RAW_H +#ifdef __cplusplus +extern "C" { +#endif + #include "lasp_config.h" #include "lasp_assert.h" #include "lasp_tracer.h" @@ -23,40 +27,6 @@ #error "LASP_USE_BLAS should be set to either 0 or 1" #endif -#if LASP_DOUBLE_PRECISION == 1 -#define c_real creal -#define c_imag cimag -#define d_abs fabs -#define c_abs cabs -#define c_conj conj -#define d_atan2 atan2 -#define d_acos acos -#define d_sqrt sqrt -#define c_exp cexp -#define d_sin sin -#define d_cos cos -#define d_pow pow -#define d_log10 log10 -#define d_ln log -#define d_epsilon (DBL_EPSILON) - -#else // LASP_DOUBLE_PRECISION not defined -#define c_conj conjf -#define c_real crealf -#define c_imag cimagf -#define d_abs fabsf -#define c_abs cabsf -#define d_atan2 atan2f -#define d_acos acosf -#define d_sqrt sqrtf -#define c_exp cexpf -#define d_sin sinf -#define d_cos cosf -#define d_pow powf -#define d_log10 log10f -#define d_ln logf -#define d_epsilon (FLT_EPSILON) -#endif // LASP_DOUBLE_PRECISION /// Positive infinite #define d_inf (INFINITY) @@ -502,6 +472,9 @@ static inline void c_conj_inplace(c res[],us size) { } #endif // LASP_USE_BLAS } +#ifdef __cplusplus +} +#endif #endif // LASP_MATH_RAW_H ////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_nprocs.h b/lasp/c/lasp_nprocs.h index 7ba5458..ea2f34c 100644 --- a/lasp/c/lasp_nprocs.h +++ b/lasp/c/lasp_nprocs.h @@ -3,16 +3,23 @@ // Author: J.A. de Jong - ASCEE // // Description: Implemententation of a function to determine the number -// of processors. +// of processors. Useful for spreading out computational load over multiple +// cores. ////////////////////////////////////////////////////////////////////// #pragma once #ifndef LASP_NPROCS_H #define LASP_NPROCS_H #include "lasp_types.h" +#ifdef __cplusplus +extern C { +#endif -/** - * @return The number of SMP processors - */ -us getNumberOfProcs(); + /** + * @return The number of SMP processors + */ + us getNumberOfProcs(); -#endif // LASP_NPROCS_H \ No newline at end of file +#ifdef __cplusplus +} +#endif +#endif // LASP_NPROCS_H diff --git a/lasp/c/lasp_siggen.c b/lasp/c/lasp_siggen.c deleted file mode 100644 index 702a527..0000000 --- a/lasp/c/lasp_siggen.c +++ /dev/null @@ -1,472 +0,0 @@ -// lasp_siggen.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// Signal generator implementation -////////////////////////////////////////////////////////////////////// -/* #define TRACERPLUS (-5) */ -#include "lasp_siggen.h" -#include "lasp_alloc.h" -#include "lasp_assert.h" -#include "lasp_mat.h" - -/** The fixed number of Newton iterations t.b.d. for tuning the sweep start and - * stop frequency in logarithmic sweeps */ -#define NITER_NEWTON 20 - -/** The number of Bytes of space for the signal-specific data in the Siggen - * structure */ -#define PRIVATE_SIZE 64 - -typedef enum { - SINEWAVE = 0, - NOISE, - SWEEP, -} SignalType; - -typedef struct Siggen { - SignalType signaltype; - d fs; // Sampling frequency [Hz] - d level_amp; - char private_data[PRIVATE_SIZE]; -} Siggen; - -typedef struct { - d curtime; - d omg; -} SinewaveSpecific; - -typedef struct { - us N; - vd data; - us index; -} PeriodicSpecific; - -typedef struct { - d V1, V2, S; - int phase; - Sosfilterbank* colorfilter; -} NoiseSpecific; - -static d level_amp(d level_dB){ - return pow(10, level_dB/20); -} - -Siggen* Siggen_create(SignalType type, const d fs,const d level_dB) { - - fsTRACE(15); - - Siggen* siggen = a_malloc(sizeof(Siggen)); - siggen->signaltype = type; - siggen->fs = fs; - siggen->level_amp = level_amp(level_dB); - - feTRACE(15); - return siggen; -} - -Siggen* Siggen_Sinewave_create(const d fs, const d freq,const d level_dB) { - fsTRACE(15); - - Siggen* sine = Siggen_create(SINEWAVE, fs, level_dB); - dbgassert(sizeof(SinewaveSpecific) <= sizeof(sine->private_data), - "Allocated memory too small"); - SinewaveSpecific* sp = (SinewaveSpecific*) sine->private_data; - sp->curtime = 0; - sp->omg = 2*number_pi*freq; - - feTRACE(15); - return sine; -} - -Siggen* Siggen_Noise_create(const d fs, const d level_dB, Sosfilterbank* colorfilter) { - fsTRACE(15); - - Siggen* noise = Siggen_create(NOISE, fs, level_dB); - dbgassert(sizeof(NoiseSpecific) <= sizeof(noise->private_data), - "Allocated memory too small"); - NoiseSpecific* wn = (NoiseSpecific*) noise->private_data; - wn->phase = 0; - wn->V1 = 0; - wn->V2 = 0; - wn->S = 0; - wn->colorfilter = colorfilter; - - feTRACE(15); - return noise; -} - -Siggen* Siggen_Sweep_create(const d fs,const d fl_,const d fu_, - const d Ts,const d Tq, const us flags, const d level_dB) { - fsTRACE(15); - - Siggen* sweep = Siggen_create(SWEEP, fs, level_dB); - dbgassert(sizeof(PeriodicSpecific) <= sizeof(sweep->private_data), - "Allocated memory too small"); - - bool forward_sweep = flags & SWEEP_FLAG_FORWARD; - bool backward_sweep = flags & SWEEP_FLAG_BACKWARD; - - // Set pointer to inplace data storage - dbgassert(!(forward_sweep && backward_sweep), "Both forward and backward flag set"); - - PeriodicSpecific* sp = (PeriodicSpecific*) sweep->private_data; - if(fl_ < 0 || fu_ < 0 || Ts <= 0) { - return NULL; - } - - - const d Dt = 1/fs; // Deltat - - // Estimate N, the number of samples in the sweep part (non-quiescent part): - const us N = (us) (Ts*fs); - const us Nq = (us) (Tq*fs); - iVARTRACE(15, N); - sp->data = vd_alloc(N+Nq); - vd* data = &(sp->data); - /* Set the last part, the quiescent tail to zero */ - dmat_set(data,0.0); - - sp->N = N+Nq; - sp->index = 0; - - // Obtain flags and expand - d phase = 0; - d fl, fu; - - /* Swap fl and fu for a backward sweep */ - if(backward_sweep) { - fu = fl_; - fl = fu_; - } - else { - /* Case of continuous sweep, or forward sweep */ - fl = fl_; - fu = fu_; - } - - /* Linear sweep */ - if(flags & SWEEP_FLAG_LINEAR) { - TRACE(15, "linear sweep"); - - if(forward_sweep || backward_sweep) { - /* Forward or backward sweep */ - TRACE(15, "Forward or backward sweep"); - us K = (us) (Dt*(fl*N+0.5*(N-1)*(fu-fl))); - d eps_num = ((d) K)/Dt - fl*N-0.5*(N-1)*(fu-fl); - d eps = eps_num/(0.5*(N-1)); - iVARTRACE(15, K); - dVARTRACE(15, eps); - - for(us n = 0; n= 0, "BUG"); - - phase += 2*number_pi*Dt*fn; - /* dVARTRACE(17, phase); */ - /* setvecval(data, n, fn); */ - /* setvecval(data, n, phase); */ - } - /* This should be a very small number!! */ - dVARTRACE(15, phase); - } - - } - else if(flags & SWEEP_FLAG_EXPONENTIAL) { - - TRACE(15, "exponential sweep"); - if(forward_sweep || backward_sweep) { - /* Forward or backward sweep */ - TRACE(15, "Forward or backward sweep"); - d k1 = (fu/fl); - us K = (us) (Dt*fl*(k1-1)/(d_pow(k1,1.0/N)-1)); - d k = k1; - - /* Iterate k to the right solution */ - d E; - for(us iter=0;iter< 10; iter++) { - E = 1 + K/(Dt*fl)*(d_pow(k,1.0/N)-1) - k; - d dEdk = K/(Dt*fl)*d_pow(k,1.0/N)/(N*k)-1; - k -= E/dEdk; - } - - iVARTRACE(15, K); - dVARTRACE(15, k1); - dVARTRACE(15, k); - dVARTRACE(15, E); - - for(us n = 0; n= 0, "BUG"); - - phase += 2*number_pi*Dt*fn; - while(phase > 2*number_pi) phase -= 2*number_pi; - /* dVARTRACE(17, phase); */ - /* setvecval(data, n, fn); */ - /* setvecval(data, n, phase); */ - } - /* This should be a very small number!! */ - dVARTRACE(15, phase); - - } - } - - feTRACE(15); - return sweep; -} - -static void Siggen_periodic_free(PeriodicSpecific* ps) { - assertvalidptr(ps); - fsTRACE(15); - vd_free(&(ps->data)); - feTRACE(15); -} -us Siggen_getN(const Siggen* siggen) { - fsTRACE(15); - assertvalidptr(siggen); - - switch(siggen->signaltype) { - case SINEWAVE: - break; - case NOISE: - break; - case SWEEP: - return ((PeriodicSpecific*) siggen->private_data)->N; - break; - default: - dbgassert(false, "Not implementend signal type"); - - } - - feTRACE(15); - return 0; -} -void Siggen_setLevel(Siggen* siggen, const d new_level_dB) { - fsTRACE(15); - - siggen->level_amp = d_pow(10, new_level_dB/20); - - feTRACE(15); -} - -void Siggen_free(Siggen* siggen) { - fsTRACE(15); - assertvalidptr(siggen); - NoiseSpecific* sp; - - switch(siggen->signaltype) { - case SWEEP: - /* Sweep specific stuff here */ - Siggen_periodic_free((PeriodicSpecific*) siggen->private_data); - break; - case SINEWAVE: - /* Sweep specific stuff here */ - break; - case NOISE: - sp = (NoiseSpecific*) siggen->private_data; - if(sp->colorfilter) { - Sosfilterbank_free(sp->colorfilter); - } - - } - - a_free(siggen); - feTRACE(15); -} - -static void Sinewave_genSignal(Siggen* siggen, SinewaveSpecific* sine, vd* samples) { - fsTRACE(10); - assertvalidptr(sine); - d ts = 1/siggen->fs; - d omg = sine->omg; - - d curtime = sine->curtime; - for(us i =0; i< samples->n_rows; i++) { - setvecval(samples, i, siggen->level_amp*sin(omg*curtime)); - curtime = curtime + ts; - } - sine->curtime = curtime; - feTRACE(10); -} - -static void Periodic_genSignal(Siggen* siggen, PeriodicSpecific* sweep, vd* samples) { - fsTRACE(10); - - for(us i=0; in_rows; i++) { - d* data = getvdval(&(sweep->data), sweep->index); - setvecval(samples, i, siggen->level_amp*(*data)); - sweep->index++; - sweep->index %= sweep->N; - } - - feTRACE(10); -} - - -static void noise_genSignal(Siggen* siggen, NoiseSpecific* wn, vd* samples) { - fsTRACE(10); - d X; - d S = wn->S; - d V1 = wn->V1; - d V2 = wn->V2; - - int phase = wn->phase; - - for(us i =0; i< samples->n_rows; i++) { - - if(wn->phase == 0) { - do { - d U1 = (d)rand() / RAND_MAX; - d U2 = (d)rand() / RAND_MAX; - - V1 = 2 * U1 - 1; - V2 = 2 * U2 - 1; - S = V1 * V1 + V2 * V2; - } while(S >= 1 || S == 0); - - X = V1 * sqrt(-2 * d_ln(S) / S); - } else - X = V2 * sqrt(-2 * d_ln(S) / S); - - phase = 1 - phase; - - setvecval(samples, i, siggen->level_amp*X); - } - if(wn->colorfilter){ - vd filtered = Sosfilterbank_filter(wn->colorfilter, - samples); - dmat_copy(samples, &filtered); - vd_free(&filtered); - } - wn->S = S; - wn->V1 = V1; - wn->V2 = V2; - wn->phase = phase; - feTRACE(10); -} - -void Siggen_genSignal(Siggen* siggen,vd* samples) { - - fsTRACE(10); - assertvalidptr(siggen); - assert_vx(samples); - - switch(siggen->signaltype) { - case SINEWAVE: - Sinewave_genSignal(siggen, - (SinewaveSpecific*) siggen->private_data, - samples); - - break; - case NOISE: - noise_genSignal(siggen, - (NoiseSpecific*) siggen->private_data, - samples); - break; - case SWEEP: - Periodic_genSignal(siggen, - (PeriodicSpecific*) siggen->private_data, - samples); - break; - default: - dbgassert(false, "Not implementend signal type"); - - } - - feTRACE(10); -} - - - -////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_siggen.h b/lasp/c/lasp_siggen.h deleted file mode 100644 index 5be6cee..0000000 --- a/lasp/c/lasp_siggen.h +++ /dev/null @@ -1,100 +0,0 @@ -// lasp_siggen.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// Header file for signal generation routines -// -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_SIGGEN_H -#define LASP_SIGGEN_H -#include "lasp_mat.h" -#include "lasp_sosfilterbank.h" - -// Define this flag to repeat a forward sweep only, or backward only. If not -// set, we do a continuous sweep -#define SWEEP_FLAG_FORWARD 1 -#define SWEEP_FLAG_BACKWARD 2 - -// Types of sweeps -#define SWEEP_FLAG_LINEAR 4 -#define SWEEP_FLAG_EXPONENTIAL 8 - -typedef struct Siggen Siggen; - -/** - * Create a sine wave signal generator - * - * @param[in] fs: Sampling frequency [Hz] - * @param[in] level: Relative level in [dB], should be between -inf and 0 - * @param[freq] Sine wave frequency [Hz] - */ -Siggen* Siggen_Sinewave_create(const d fs,const d freq,const d level_dB); - -/** - * Create a Noise signal generator. If no Sosfilterbank is provided, it will - * create white noise. Otherwise, the noise is 'colored' using the filterbank - * given in the constructor. Note that the pointer to this filterbank is - * *STOLEN*!. - * - * @param[in] fs: Sampling frequency [Hz] - * @param[in] level_dB: Relative level [dB] - * @param[in] - * - * @return Siggen* handle - */ -Siggen* Siggen_Noise_create(const d fs, const d level_dB, Sosfilterbank* colorfilter); - -/** - * Set the level of the signal generator - * @param[in] Siggen* Signal generator handle - * - * @param[in] new_level_dB The new level, in dBFS - */ -void Siggen_setLevel(Siggen*, const d new_level_dB); - - -/** - * Obtain the repetition period for a periodic excitation. - * @param[in] Siggen* Signal generator handle - * - * @param[out] N The amount of samples in one period, returns 0 if the signal - * does not repeat. - */ -us Siggen_getN(const Siggen*); - -/** - * Create a forward sweep - * - * @param[in] fs: Sampling frequency [Hz] - * @param[in] fl: Lower frequency [Hz] - * @param[in] fl: Upper frequency [Hz] - * @param[in] Ts: Sweep time [s] - * @param[in] Tq: Quescent tail time [s]. Choose this value long enough to - * avoid temporal aliasing in case of measuring impulse responses. - * @param[in] sweep_flags: Sweep period [s] - * @param[in] level: Relative level in [dB], should be between -inf and 0 - * @return Siggen* handle - */ -Siggen* Siggen_Sweep_create(const d fs,const d fl,const d fu, - const d Ts, const d Tq, const us sweep_flags, - const d level); - -/** - * Obtain a new piece of signal - * - * @param[in] Siggen* Signal generator handle - * @param[out] samples Samples to fill. Vector should be pre-allocated, but - * values will be overwritten. - */ -void Siggen_genSignal(Siggen*,vd* samples); -/** - * Free Siggen data - * - * @param[in] Siggen* Signal generator private data - */ -void Siggen_free(Siggen*); - -#endif //LASP_SIGGEN_H -////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_tracer.c b/lasp/c/lasp_tracer.c deleted file mode 100644 index d404776..0000000 --- a/lasp/c/lasp_tracer.c +++ /dev/null @@ -1,78 +0,0 @@ -// lasp_tracer.c -// -// last-edit-by: J.A. de Jong -// -// Description: -// Debug tracing code implementation -////////////////////////////////////////////////////////////////////// -#include "lasp_config.h" -#if TRACER == 1 -#include -#include "lasp_tracer.h" -#include "lasp_types.h" - -#ifdef _REENTRANT - -static __thread us ASCEE_FN_LEVEL = 0; - -static int TRACERNAME = DEFAULTTRACERLEVEL; - - -void setTracerLevel(int level) { - __sync_lock_test_and_set(&TRACERNAME, level); -} -int getTracerLevel() { - return __sync_fetch_and_add(&TRACERNAME, 0); -} - -#else - -int TRACERNAME; -static us ASCEE_FN_LEVEL = 0; - -/* setTracerLevel and getTracerLevel are defined as macros in - * tracer.h */ -#endif - -void indent_trace() { - for(us i=0;i -#include -#include -#ifdef __cplusplus -extern "C" { -#endif - -static inline void clearScreen() { - printf("\033c\n"); -} - -// Some console colors -#define RESET "\033[0m" -#define BLACK "\033[30m" /* Black */ -#define RED "\033[31m" /* Red */ -#define GREEN "\033[32m" /* Green */ -#define YELLOW "\033[33m" /* Yellow */ -#define BLUE "\033[34m" /* Blue */ -#define MAGENTA "\033[35m" /* Magenta */ -#define CYAN "\033[36m" /* Cyan */ -#define WHITE "\033[37m" /* White */ -#define BOLDBLACK "\033[1m\033[30m" /* Bold Black */ -#define BOLDRED "\033[1m\033[31m" /* Bold Red */ -#define BOLDGREEN "\033[1m\033[32m" /* Bold Green */ -#define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */ -#define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */ -#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */ -#define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */ -#define BOLDWHITE "\033[1m\033[37m" /* Bold White */ - -// Not so interesting part -#define rawstr(x) #x -#define namestr(x) rawstr(x) -#define annestr(x) namestr(x) -#define FILEWITHOUTPATH ( strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : strrchr(__FILE__, '\\') ? strrchr(__FILE__, '\\') + 1 : __FILE__ ) -// #define POS annestr(FILEWITHOUTPATH) ":" # __LINE__ << ": " -// End not so interesting part - -/** - * Produce a debug warning - */ -#define DBGWARN(a) \ - printf(RED); \ - printf("%s(%d): ", \ - __FILE__, \ - __LINE__ \ - ); \ - printf(a); \ - printf(RESET "\n"); - -/** - * Produce a runtime warning - */ -#define WARN(a) \ - printf(RED); \ - printf("WARNING: "); \ - printf(a); \ - printf(RESET "\n"); - -/** - * Fatal error, abort execution - */ -#define FATAL(a) \ - WARN(a); \ - abort(); - - - - -// **************************************** Tracer code -#ifndef TRACERPLUS -#define TRACERPLUS (0) -#endif - -// If PP variable TRACER is not defined, we automatically set it on. -#ifndef TRACER -#define TRACER 1 -#endif - - -#if TRACER == 1 -#ifndef TRACERNAME - -#ifdef __GNUC__ -#warning TRACERNAME name not set, sol TRACERNAME set to 'defaulttracer' -#else -#pragma message("TRACERNAME name not set, sol TRACERNAME set to defaulttracer") -#endif - -#define TRACERNAME defaulttracer -#endif // ifndef TRACERNAME - -/** - * Indent the rule for tracing visibility. - */ -void indent_trace(); - -// Define this preprocessor definition to overwrite -// Use -O flag for compiler to remove the dead functions! -// In that case all cout's for TRACE() are removed from code -#ifndef DEFAULTTRACERLEVEL -#define DEFAULTTRACERLEVEL (15) -#endif - -#ifdef _REENTRANT -/** - * Set the tracer level at runtime - * - * @param level - */ -void setTracerLevel(int level); - -/** - * Obtain the tracer level - * - * @return level - */ -int getTracerLevel(); - -#else // Not reentrant -extern int TRACERNAME; -#define setTracerLevel(a) TRACERNAME = a; -static inline int getTracerLevel() { return TRACERNAME;} -#endif - -#include "lasp_types.h" - -// Use this preprocessor command to introduce one TRACERNAME integer per unit -/* Introduce one static logger */ -// We trust that the compiler will eliminate 'dead code', which means -// that if variable BUILDINTRACERLEVEL is set, the inner if statement -// will not be reached. -void trace_impl(const char* pos,int line,const char * string); -void fstrace_impl(const char* file,int pos,const char* fn); -void fetrace_impl(const char* file,int pos,const char* fn); -void dvartrace_impl(const char* pos,int line,const char* varname,d var); -void cvartrace_impl(const char* pos,int line,const char* varname,c var); -void ivartrace_impl(const char* pos,int line,const char* varname,int var); -void uvartrace_impl(const char* pos,int line,const char* varname,size_t var); - -/** - * Print a trace string - */ -#define TRACE(level,trace_string) \ - if (level+TRACERPLUS>=getTracerLevel()) \ - { \ - trace_impl(FILEWITHOUTPATH,__LINE__,trace_string ); \ - } - -#define SFSG TRACE(100,"SFSG") - -/** - * Print start of function string - */ -#define fsTRACE(level) \ - if (level+TRACERPLUS>=getTracerLevel()) \ - { \ - fstrace_impl(FILEWITHOUTPATH,__LINE__, __FUNCTION__ ); \ - } - -/** - * Print end of function string - */ -#define feTRACE(level) \ - if (level+TRACERPLUS>=getTracerLevel()) \ - { \ - fetrace_impl(FILEWITHOUTPATH,__LINE__, __FUNCTION__ ); \ - } - - -/** - * Trace an int variable - */ -#define iVARTRACE(level,trace_var) \ - if (level+TRACERPLUS>=getTracerLevel()) \ - { \ - ivartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ - } - -/** - * Trace an unsigned int variable - */ -#define uVARTRACE(level,trace_var) \ - if (level+TRACERPLUS>=getTracerLevel()) \ - { \ - uvartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ - } -/** - * Trace a floating point value - */ -#define dVARTRACE(level,trace_var) \ - if (level+TRACERPLUS>=getTracerLevel()) \ - { \ - dvartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ - } -/** - * Trace a complex floating point value - */ -#define cVARTRACE(level,trace_var) \ - if (level+TRACERPLUS>=getTracerLevel()) \ - { \ - cvartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ - } - - -#else // TRACER !=1 -#define TRACE(l,a) -#define fsTRACE(l) -#define feTRACE(l) -#define setTracerLevel(a) -#define getTracerLevel() - -#define iVARTRACE(level,trace_var) -#define uVARTRACE(level,trace_var) -#define dVARTRACE(level,trace_var) -#define cVARTRACE(level,trace_var) - - -#endif // ######################################## TRACER ==1 - - -#ifdef __cplusplus -} -#endif -#endif // LASP_TRACER_H -////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_window.c b/lasp/c/lasp_window.c deleted file mode 100644 index 9f11123..0000000 --- a/lasp/c/lasp_window.c +++ /dev/null @@ -1,111 +0,0 @@ -// window.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// - -#include "lasp_window.h" -#include "lasp_signals.h" -#include -/** - * Compute the Hann window - * - * @param i index - * @param N Number of indices - */ -static d hann(us i,us N) { - dbgassert(in_rows; - d (*win_fun)(us,us); - switch (wintype) { - case Hann: { - win_fun = hann; - break; - } - case Hamming: { - win_fun = hamming; - break; - } - case Rectangular: { - win_fun = rectangle; - break; - } - case Bartlett: { - win_fun = bartlett; - break; - } - case Blackman: { - win_fun = blackman; - break; - } - default: - DBGWARN("BUG: Unknown window function"); - abort(); - break; - } - us index; - for(index=0;index #include "lasp_config.h" #include "lasp_daqdata.h" #include "lasp_deviceinfo.h" #include "lasp_types.h" #include #include +#include + +/** + * @brief Information regarding a stream. + */ +class StreamStatus { +public: + bool isRunning = false; + bool error = false; + std::string error_msg{}; + + /** + * @brief Returns true if everything is OK with a certain stream and the + * stream is running. + * + * @return as described above. + */ + bool runningOK() const { return isRunning && !error; } +}; /** * @brief Callback of DAQ for input data. Callback should return * false for a stop request. */ -using InDaqCallback = std::function; +using InDaqCallback = std::function; /** - * @brief + * @brief */ -using OutDaqCallback = std::function; +using OutDaqCallback = std::function; /** * @brief Base cass for all DAQ instances @@ -48,7 +66,7 @@ public: * data is presented, the function is called with a nullptr as argument. */ virtual void start(std::optional inCallback, - std::optional outCallback); + std::optional outCallback) = 0; /** * @brief Stop the Daq device. Throws an exception if the device is not @@ -63,7 +81,7 @@ public: */ virtual bool isRunning() const = 0; - virtual ~Daq(){}; + virtual ~Daq() = default; /** * @brief Returns the number of enabled input channels * @@ -118,6 +136,14 @@ public: us framesPerBlock() const { return availableFramesPerBlock.at(framesPerBlockIndex); } + + /** + * @brief Get stream status corresponding to current DAQ. + * + * @return StreamStatus object. + */ + virtual StreamStatus getStreamStatus() const = 0; + /** * @brief Whether the device runs in duplex mode (both input and output), or * false if only input / only output. diff --git a/lasp/device/lasp_daqconfig.h b/lasp/device/lasp_daqconfig.h index 5b233a1..7b161cc 100644 --- a/lasp/device/lasp_daqconfig.h +++ b/lasp/device/lasp_daqconfig.h @@ -1,8 +1,8 @@ #pragma once -#include "lasp_config.h" -#include "lasp_types.h" #include #include +#include "lasp_config.h" +#include "lasp_types.h" using std::string; using std::to_string; @@ -200,7 +200,7 @@ public: * @param deviceinfo DeviceInfo structure */ DaqConfiguration(const DeviceInfo &DeviceInfo); - DaqConfiguration() = delete; + DaqConfiguration(); /** * @brief Check to see whether the DAQ configuration matches with the device. diff --git a/lasp/device/lasp_daqdata.h b/lasp/device/lasp_daqdata.h index fbbf94c..fd71441 100644 --- a/lasp/device/lasp_daqdata.h +++ b/lasp/device/lasp_daqdata.h @@ -1,10 +1,10 @@ #pragma once -#include "lasp_daqconfig.h" -#include "lasp_types.h" #include #include #include #include +#include "lasp_daqconfig.h" +#include "lasp_types.h" /** * @brief Data coming from / going to DAQ. **Non-interleaved format**, which diff --git a/lasp/device/lasp_devicepybind.cpp b/lasp/device/lasp_devicepybind.cpp index 481c318..71fd668 100644 --- a/lasp/device/lasp_devicepybind.cpp +++ b/lasp/device/lasp_devicepybind.cpp @@ -1,31 +1,33 @@ -#if 0 -#include "lasp_devicepybind.h" #include #include -#include #include #include -#include +#include +#include "lasp_daqconfig.h" +#include "lasp_deviceinfo.h" using std::cerr; namespace py = pybind11; -PYBIND11_MODULE(lasp_daq, m) { +PYBIND11_MODULE(lasp_device, m) { - m.doc() = "Lasp DAQ interface"; + m.doc() = "Lasp device interface"; /// DataType - py::class_ datatype(m, "DataType"); - datatype.def_readonly("name", &DataType::name); - datatype.def_readonly("sw", &DataType::sw); - datatype.def_readonly("is_floating", &DataType::is_floating); + py::class_ dtype_desc(m, "DataTypeDescriptor"); + dtype_desc.def_readonly("name", &DataTypeDescriptor::name); + dtype_desc.def_readonly("sw", &DataTypeDescriptor::sw); + dtype_desc.def_readonly("is_floating", &DataTypeDescriptor::is_floating); + + py::enum_(dtype_desc, "DataType") + .value("dtype_fl32", DataTypeDescriptor::DataType::dtype_fl32) + .value("dtype_fl64", DataTypeDescriptor::DataType::dtype_fl64) + .value("dtype_int8", DataTypeDescriptor::DataType::dtype_int8) + .value("dtype_int16", DataTypeDescriptor::DataType::dtype_int16) + .value("dtype_int32", DataTypeDescriptor::DataType::dtype_int32).export_values(); + + dtype_desc.def_readonly("dtype", &DataTypeDescriptor::dtype); - m.attr("dtype_fl32") = dtype_fl32; - m.attr("dtype_fl64") = dtype_fl64; - m.attr("dtype_int8") = dtype_int8; - m.attr("dtype_int16") = dtype_int16; - m.attr("dtype_int24") = dtype_int24; - m.attr("dtype_int32") = dtype_int32; /// DaqApi py::class_ daqapi(m, "DaqApi"); @@ -105,132 +107,5 @@ PYBIND11_MODULE(lasp_daq, m) { daqconfig.def("match", &DaqConfiguration::match); - /// Daq - py::class_ daq(m, "Daq"); - daq.def(py::init()); - /* daq.def("start */ - /* daqconfig.def(py::init< */ } -const d QUEUE_BUFFER_TIME = 0.1; - -PyDaq::PyDaq(const DeviceInfo &devinfo, const DaqConfiguration &config) { - - _daq = Daq::createDaq(devinfo, config); - - _stopThread = false; - _threadReady = false; -} -PyDaq::~PyDaq() { - if (_daqThread) { - } -} - -/* py::array arrayFromBuffer(void *buf, us n_rows, us n_cols, DataType dtype) { */ -/* // https://github.com/pybind/pybind11/issues/27 */ -/* py::size_t dims[] = {n_rows, n_cols}; */ -/* py::buffer_info; */ -/* switch (dtype) { */ -/* case dtype_int8: */ -/* break; */ -/* case dtype_int16: */ -/* break; */ -/* case dtype_int24: */ -/* break; */ -/* case dtype_fl32: */ -/* break; */ -/* case dtype_fl64: */ -/* break; */ -/* } */ -/* /1* buffer_info = py::buffer_info( *1/ */ -/* /1* py::buffer_info(buf, sizeof(float), *1/ */ -/* /1* py::format_descriptor::format(), *1/ */ -/* /1* 2, *1/ */ -/* /1* dims, (float*) buf));; *1/ */ -/* /1* } *1/ */ -/* py::array res; */ -/* } */ - -void PyDaq::start(py::function audioCallback) { - - // Number of input channels, including monitor channel - const us neninchannels = _daq->neninchannels(true); - const us nenoutchannels = _daq->nenoutchannels(); - const bool input = neninchannels > 0; - const bool output = nenoutchannels > 0; - - SafeQueue *inQueue = input ? &_inQueue : nullptr; - SafeQueue *outQueue = output ? &_outQueue : nullptr; - - _daq->start(inQueue, outQueue); - - _audioCallback = audioCallback; - - _daqThread = new std::thread(&PyDaq::threadFcn, this); - - // Wait unitil thread is ready before returning - do { - std::this_thread::sleep_for(std::chrono::milliseconds(1)); - } while (!_threadReady); -} - -void PyDaq::stop() { - _daq->stop(); - _stopThread = true; - _threadReady = false; - _daqThread->join(); - delete _daqThread; - _daqThread = nullptr; -} - -void PyDaq::threadFcn() { - - assert(_daq); - - void *inbuffer = nullptr; - void *outbuffer = nullptr; - - const double sleeptime = _daq->framesPerBlock() / (8 * _daq->samplerate()); - // Sleep time in microseconds - const us sleeptime_us = (us)(sleeptime * 1e6); - - const us nblocks_buffer = (us)std::max( - 1, QUEUE_BUFFER_TIME * _daq->samplerate() / _daq->framesPerBlock()); - - auto descr = _daq->dtypeDescr(); - const us nBytesPerChan = _daq->framesPerBlock() * descr.sw; - - _threadReady = true; - - // Number of input channels, including monitor channel - const us neninchannels = _daq->neninchannels(true); - const us nenoutchannels = _daq->nenoutchannels(); - const bool input = neninchannels > 0; - const bool output = nenoutchannels > 0; - - cerr << "Check if it is required to call import_array()" << endl; - /* { */ - /* py::gil_scoped_acquire aq; */ - /* import_array(void); */ - - /* } */ - _threadReady = true; - if (output) { - for (us i = 0; i < nblocks_buffer; i++) { - outbuffer = new uint8_t[nBytesPerChan * nenoutchannels]; - memset(outbuffer, 0, nBytesPerChan * nenoutchannels); - _outQueue.enqueue(outbuffer); - } - outbuffer = nullptr; - } - - while (!_stopThread) { - if (output) { - py::gil_scoped_acquire gil; - while (_outQueue.size() < nblocks_buffer) { - /* py::array( */ - } - } - } // End of while loop, stopthread -} -#endif diff --git a/lasp/device/lasp_devicepybind.h b/lasp/device/lasp_devicepybind.h deleted file mode 100644 index ef4a790..0000000 --- a/lasp/device/lasp_devicepybind.h +++ /dev/null @@ -1,44 +0,0 @@ -#pragma once -#include "lasp_daq.h" -#include -#include - -#include -namespace py = pybind11; - -namespace std { -class thread; -}; - -class PyDaq { - Daq *_daq = nullptr; - std::thread *_daqThread = nullptr; - - SafeQueue _inQueue; - SafeQueue _outQueue; - - py::function _audioCallback; - - std::atomic _stopThread; - std::atomic _threadReady; - -public: - PyDaq(const DeviceInfo &devinfo, const DaqConfiguration &config); - ~PyDaq(); - /** - * @brief Start the Daq data acquisition - * @param audioCallback Python callback which is called with in and/or out - data of the daq - * as arguments - - */ - void start(py::function audioCallback); - void stop(); - bool isRunning() { return _daq->isRunning(); } - - static std::vector getDeviceInfo(); - const Daq &daq() const { return *_daq; } - -private: - void threadFcn(); -}; diff --git a/lasp/device/lasp_rtaudiodaq.cpp b/lasp/device/lasp_rtaudiodaq.cpp index e0427f3..c60fff3 100644 --- a/lasp/device/lasp_rtaudiodaq.cpp +++ b/lasp/device/lasp_rtaudiodaq.cpp @@ -300,14 +300,10 @@ public: return 0; } - ~RtAudioDaq() { - if (isRunning()) { - stop(); - } - if (rtaudio.isStreamOpen()) { - rtaudio.closeStream(); - } - } + // RtAudio documentation says: if a stream is open, it will be stopped and + // closed automatically on deletion. Therefore the destructor here is a + // default one. + ~RtAudioDaq() = default; }; std::unique_ptr createRtAudioDevice(const DeviceInfo &devinfo, diff --git a/lasp/device/lasp_streammgr.cpp b/lasp/device/lasp_streammgr.cpp new file mode 100644 index 0000000..656a3b7 --- /dev/null +++ b/lasp/device/lasp_streammgr.cpp @@ -0,0 +1,151 @@ +#include +#include +#include "lasp_streammgr.h" + +InDataHandler::InDataHandler(StreamMgr &mgr) : _mgr(mgr) { + mgr.addInDataHandler(*this); +} +InDataHandler::~InDataHandler() { _mgr.removeInDataHandler(*this); } + + +StreamMgr &StreamMgr::getInstance() { + static StreamMgr mgr; + return mgr; +} + +bool StreamMgr::inCallback(const DaqData &data) { + std::scoped_lock lck(_inDataHandler_mtx); + + for (auto &handler : _inDataHandlers) { + bool res = handler->inCallback(data); + if (!res) + return false; + } + return true; +} + +/** + * @brief Converts from double precision floating point to output signal in + * non-interleaving format. + * + * @tparam T + * @param data + * @param signal + * + * @return + */ +template bool fillData(DaqData &data, const vd &signal) { + assert(data.nframes == signal.size()); + + T* res = reinterpret_cast(data.raw_ptr()); + if (std::is_floating_point()) { + for (us ch = 0; ch < data.nchannels; ch++) { + for (us frame = 0; frame < data.nframes; frame++) { + res[ch * data.nframes + frame ] = signal[frame]; + } + } + } else { + for (us ch = 0; ch < data.nchannels; ch++) { + for (us frame = 0; frame < data.nframes; frame++) { + res[ch * data.nframes + frame] = + (signal[frame] * std::numeric_limits::max()); + } + } + } + + return true; +} + +void StreamMgr::setSiggen(std::shared_ptr siggen) { + + std::scoped_lock lck(_siggen_mtx); + _siggen = siggen; + // If not set to nullptr, and a stream is running, we update the signal + // generator by resetting it. + if(isStreamRunning(StreamType::outputType) && siggen) { + const Daq* daq = getDaq(StreamType::outputType); + assert(daq != nullptr); + // Reset the signal generator. + _siggen->reset(daq->samplerate()); + } + +} + +bool StreamMgr::outCallback(DaqData& data) { + std::scoped_lock lck(_siggen_mtx); + if(_siggen) { + vd signal = _siggen->genSignal(data.nframes); + switch (data.dtype) { + case (DataTypeDescriptor::DataType::dtype_fl32): + fillData(data, signal); + break; + case (DataTypeDescriptor::DataType::dtype_fl64): + fillData(data, signal); + break; + case (DataTypeDescriptor::DataType::dtype_int8): + fillData(data, signal); + break; + case (DataTypeDescriptor::DataType::dtype_int16): + fillData(data, signal); + break; + case (DataTypeDescriptor::DataType::dtype_int32): + fillData(data, signal); + break; + } + } else { + // Set all values to 0. + std::fill(data.raw_ptr(), data.raw_ptr()+data.size_bytes(), 0); + } + return true; +} + +StreamMgr::~StreamMgr() { stopAllStreams(); } +void StreamMgr::stopAllStreams() { _streams.clear(); } + +void StreamMgr::startStream(const DeviceInfo &devinfo, + const DaqConfiguration &config) { + + bool isInput = std::count(config.eninchannels.begin(), + config.eninchannels.end(), true) > 0; + bool isOutput = std::count(config.enoutchannels.begin(), + config.enoutchannels.end(), true) > 0; + bool isDuplex = isInput && isOutput; + + if (!isInput && !isOutput) { + throw std::runtime_error( + "Neither input, nor output channels enabled for stream"); + } + StreamType streamtype; + + if (isDuplex) { + if (_streams.size()) { + throw std::runtime_error("Error: a stream is already running. Please " + "first stop existing stream"); + } + streamtype = StreamType::duplex; + + } else if (isInput) { + if (_streams[StreamType::input]) { + throw std::runtime_error("Error: input stream is already running. Please " + "first stop existing stream"); + } + streamtype = StreamType::input; + } else if (isOutput) { + if (_streams[StreamType::output]) { + throw std::runtime_error( + "Error: output stream is already running. Please " + "first stop existing stream"); + } + streamtype = StreamType::input; + } + _streams[streamtype] = Daq::createDaq(devinfo, config); +} + +void StreamMgr::addInDataHandler(InDataHandler &handler) { + std::scoped_lock lck(_inDataHandler_mtx); + _inDataHandlers.push_back(&handler); +} +void StreamMgr::removeInDataHandler(InDataHandler &handler) { + std::scoped_lock lck(_inDataHandler_mtx); + _inDataHandlers.remove(&handler); +} diff --git a/lasp/device/lasp_streammgr.h b/lasp/device/lasp_streammgr.h new file mode 100644 index 0000000..b6eec04 --- /dev/null +++ b/lasp/device/lasp_streammgr.h @@ -0,0 +1,187 @@ +#pragma once +#include "lasp_daq.h" +#include "lasp_siggen.h" +#include +#include +#include + +class StreamMgr; + +/** + * @brief Information regarding a stream. + */ +class StreamStatus { +public: + bool isRunning = false; + bool error = false; + std::string error_msg{}; + + /** + * @brief Returns true if everything is OK with a certain stream and the + * stream is running. + * + * @return as described above. + */ + bool runningOK() const { return isRunning && !error; } +}; + +class InDataHandler { + +protected: + StreamMgr &_mgr; + +public: + virtual ~InDataHandler(); + + /** + * @brief When constructed, the handler is added to the stream manager, which + * will call the handlers's inCallback() until the point of destruction. + * + * @param mgr Stream manager. + */ + InDataHandler(StreamMgr &mgr); + + /** + * @brief This function is called when input data from a DAQ is available. + * + * @param daqdata Input data from DAQ + * + * @return true if no error. False to stop the stream from running. + */ + virtual bool inCallback(const DaqData &daqdata) = 0; +}; + +/** + * @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. + */ +class StreamMgr { + + enum class StreamType : us { + /** + * @brief Input-only stream + */ + input = 1 << 0, + /** + * @brief Output-only stream + */ + output = 1 << 1, + /** + * @brief Duplex stream (does both input and output) + */ + duplex = 1 << 2, + /** + * @brief Either input, or duplex + */ + inputType = 1 << 3, + /** + * @brief Either output, or duplex + */ + outputType = 1 << 4 + }; + + /** + * @brief Storage for streams + */ + std::map> _streams; + + /** + * @brief All indata handlers are called when input data is available. Note + * that they can be called from different threads and should take care of + * thread-safety. + */ + std::list _inDataHandlers; + std::mutex _inDataHandler_mtx; + + /** + * @brief Signal generator in use to generate output data. Currently + * implemented as to generate the same data for all output channels. + */ + std::shared_ptr _siggen; + std::mutex _siggen_mtx; + + StreamMgr(); + + friend class InDataHandler; + friend class Siggen; + +public: + StreamMgr(const StreamMgr &) = delete; + StreamMgr &operator=(const StreamMgr &) = delete; + ~StreamMgr(); + + /** + * @brief Get access to stream manager instance + * + * @return Reference to stream manager. + */ + static StreamMgr &getInstance(); + + /** + * @brief Start a stream. + * + * @param devinfo Device information of device + * @param config Configuration of device + */ + void startStream(const DeviceInfo &devinfo, const DaqConfiguration &config); + + /** + * @brief Check if a certain stream is running. If running with no errors, it + * returns true. If an error occured, or the stream is not running, it gives + * false. + * + * @param type The type of stream to check for. + */ + bool isStreamRunningOK(const StreamType type) const { + return getStreamStatus(type).runningOK(); + } + + /** + * @brief Get the streamstatus object corresponding to a given stream. + * + * @param type Type of stream, input, inputType, etc. + * + * @return status object. + */ + StreamStatus getStreamStatus(const StreamType type) const; + + /** + * @brief Get DAQ pointer for a given stream. Gives a nullptr if stream is + * not running. + * + * @param type The stream type to get a DAQ ptr for. + * + * @return Pointer to DAQ + */ + const Daq *getDaq(StreamType type) const; + + /** + * @brief Stop stream of given type (input / output/ duplex); + * + * @param stype The stream type to stop. + */ + void stopStream(StreamType stype); + + /** + * @brief Stop and delete all streams. Also called on destruction of the + * StreamMgr. + */ + void stopAllStreams(); + +private: + bool inCallback(const DaqData &data); + bool outCallback(DaqData &data); + + void addInDataHandler(InDataHandler &handler); + void removeInDataHandler(InDataHandler &handler); + + /** + * @brief Set active signal generator for output streams. Only one `Siggen' + * is active at the same time. Siggen controls its own data race protection + * using a mutex. + * + * @param s Siggen pointer + */ + void setSiggen(std::shared_ptr s); +}; diff --git a/lasp/device/lasp_uldaq.cpp b/lasp/device/lasp_uldaq.cpp index dbad5bc..87c127f 100644 --- a/lasp/device/lasp_uldaq.cpp +++ b/lasp/device/lasp_uldaq.cpp @@ -8,7 +8,6 @@ #include #include #include -#include #include #include #include diff --git a/lasp/dsp/CMakeLists.txt b/lasp/dsp/CMakeLists.txt new file mode 100644 index 0000000..344f028 --- /dev/null +++ b/lasp/dsp/CMakeLists.txt @@ -0,0 +1,21 @@ +configure_file(lasp_config.h.in lasp_config.h) + +set(lasp_dsp_files + lasp_filter.cpp + lasp_siggen.cpp + lasp_siggen_impl.cpp + lasp_window.cpp + ) + +SET_SOURCE_FILES_PROPERTIES(${lasp_dsp_files} PROPERTIES COMPILE_FLAGS + "-DARMA_DONT_USE_WRAPPER") + + +add_library(lasp_dsp_lib STATIC ${lasp_dsp_files}) + +pybind11_add_module(lasp_dsp lasp_dsp_pybind.cpp) + +target_link_libraries(lasp_dsp PRIVATE lasp_dsp_lib ${BLAS_LIBRARIES}) + +target_include_directories(lasp_dsp_lib PUBLIC ${CMAKE_CURRENT_BINARY_DIR}) + diff --git a/lasp/dsp/__init__.py b/lasp/dsp/__init__.py new file mode 100644 index 0000000..47e3ef5 --- /dev/null +++ b/lasp/dsp/__init__.py @@ -0,0 +1 @@ +from .lasp_dsp import * diff --git a/lasp/c/lasp_config.h.in b/lasp/dsp/lasp_config.h.in similarity index 100% rename from lasp/c/lasp_config.h.in rename to lasp/dsp/lasp_config.h.in diff --git a/lasp/dsp/lasp_dsp_pybind.cpp b/lasp/dsp/lasp_dsp_pybind.cpp new file mode 100644 index 0000000..654f82b --- /dev/null +++ b/lasp/dsp/lasp_dsp_pybind.cpp @@ -0,0 +1,20 @@ +#include "lasp_window.h" +#include +#include + +using std::cerr; +namespace py = pybind11; + +PYBIND11_MODULE(lasp_dsp, m) { + + py::class_ w(m, "Window"); + + py::enum_(w, "WindowType") + .value("Hann", Window::WindowType::Hann) + .value("Hamming", Window::WindowType::Hamming) + .value("Bartlett", Window::WindowType::Bartlett) + .value("Blackman", Window::WindowType::Bartlett) + .value("Rectangular", Window::WindowType::Rectangular) + .export_values(); + +} diff --git a/lasp/dsp/lasp_filter.cpp b/lasp/dsp/lasp_filter.cpp new file mode 100644 index 0000000..e69de29 diff --git a/lasp/dsp/lasp_filter.h b/lasp/dsp/lasp_filter.h new file mode 100644 index 0000000..1c1c48d --- /dev/null +++ b/lasp/dsp/lasp_filter.h @@ -0,0 +1,14 @@ +#pragma once +#include +#include "lasp_types.h" +#include "lasp_mathtypes.h" +/** + * @brief Filter used to pre-filter a double-precision floating point data + * stream. + */ +class Filter { +public: + virtual void filter(const vd &inout) = 0; + virtual ~Filter() = 0; +}; + diff --git a/lasp/dsp/lasp_mathtypes.h b/lasp/dsp/lasp_mathtypes.h new file mode 100644 index 0000000..6ecce1d --- /dev/null +++ b/lasp/dsp/lasp_mathtypes.h @@ -0,0 +1,45 @@ +#pragma once +#include +#include "lasp_types.h" +#include + +#if LASP_DOUBLE_PRECISION == 1 +#define c_real creal +#define c_imag cimag +#define d_abs fabs +#define c_abs cabs +#define c_conj conj +#define d_atan2 atan2 +#define d_acos acos +#define d_sqrt sqrt +#define c_exp cexp +#define d_sin sin +#define d_cos cos +#define d_pow pow +#define d_log10 log10 +#define d_ln log +#define d_epsilon (DBL_EPSILON) + +#else // LASP_DOUBLE_PRECISION not defined +#define c_conj conjf +#define c_real crealf +#define c_imag cimagf +#define d_abs fabsf +#define c_abs cabsf +#define d_atan2 atan2f +#define d_acos acosf +#define d_sqrt sqrtf +#define c_exp cexpf +#define d_sin sinf +#define d_cos cosf +#define d_pow powf +#define d_log10 log10f +#define d_ln logf +#define d_epsilon (FLT_EPSILON) + +#endif // LASP_DOUBLE_PRECISION + +using vd = arma::Col; +using dmat = arma::Mat; + +const d number_pi = arma::datum::pi; diff --git a/lasp/dsp/lasp_siggen.cpp b/lasp/dsp/lasp_siggen.cpp new file mode 100644 index 0000000..1e74adb --- /dev/null +++ b/lasp/dsp/lasp_siggen.cpp @@ -0,0 +1,42 @@ +#include "lasp_siggen.h" +#include "lasp_mathtypes.h" +#include +#include + +inline d level_amp(d level_dB){ + return pow(10, level_dB/20); +} + +using mutexlock = std::scoped_lock; + +vd Siggen::genSignal(const us nframes) { + mutexlock lck(_mtx); + + vd signal(nframes, arma::fill::value(_dc_offset)); + if (!_muted) { + vd signal_dynamic = _level_linear*genSignalUnscaled(nframes); + if(_filter) { + _filter->filter(signal_dynamic); + } + signal += signal_dynamic; + } + + return signal; +} +void Siggen::setFilter(std::shared_ptr& filter) { + mutexlock lck(_mtx); + _filter = filter; +} +void Siggen::setDCOffset(const d offset) { + mutexlock lck(_mtx); + _dc_offset = offset; +} +void Siggen::setLevel(const d level,bool dB) { + mutexlock lck(_mtx); + _level_linear = dB ? level_amp(level) : level; +} +void Siggen::reset(const d newFs) { + mutexlock lck(_mtx); + fs = newFs; + resetImpl(); +} diff --git a/lasp/dsp/lasp_siggen.h b/lasp/dsp/lasp_siggen.h new file mode 100644 index 0000000..f95fe6c --- /dev/null +++ b/lasp/dsp/lasp_siggen.h @@ -0,0 +1,77 @@ +#pragma once +#include +#include "lasp_types.h" +#include "lasp_mathtypes.h" +#include "lasp_filter.h" + +class StreamMgr; +class DaqData; + +/** + * @brief Signal generation base class + */ +class Siggen { +private: + std::shared_ptr _filter; + d _dc_offset = 0, _level_linear = 1; + bool _muted = false; + std::mutex _mtx; + +protected: + d fs = -1; + + virtual void resetImpl() = 0; + virtual vd genSignalUnscaled(const us nframes) = 0; + +public: + virtual ~Siggen() = default; + + /** + * @brief Set a post-filter on the signal. For example to EQ the signal, or + * otherwise to shape the spectrum. + * + * @param f The filter to install. + */ + void setFilter(std::shared_ptr& f); + + /** + * @brief Set a linear DC offset value to the signal + * + * @param offset + */ + void setDCOffset(d offset); + + /** + * @brief Mute the signal. Passes through the DC offset. + * + * @param mute if tre + */ + void setMute(bool mute = true) { _muted = mute; } + + /** + * @brief Set the level of the signal generator + * + * @param level The new level. If dB == true, it is treated as a level, and + * pow(10, level/20) is installed as the linear gain. + * @param dB if false, level is treated as linear gain value. + */ + void setLevel(const d level, bool dB=true); + + /** + * @brief Reset the signal generator. Should be called whenever the output is + * based on a different sampling frequency. Note that derived classes from + * this class should call the base class! + * + * @param newFs New sampling frequency to use. + */ + void reset(const d newFs); + + /** + * @brief Called whenever the implementation needs to create new samples. + * + * @param nframes + * + * @return Array of samples with length nframes + */ + vd genSignal(const us nframes); +}; diff --git a/lasp/dsp/lasp_siggen_impl.cpp b/lasp/dsp/lasp_siggen_impl.cpp new file mode 100644 index 0000000..c76337d --- /dev/null +++ b/lasp/dsp/lasp_siggen_impl.cpp @@ -0,0 +1,245 @@ +// lasp_siggen_impl.cpp +// +// Author: J.A. de Jong -ASCEE +// +// Description: +// Signal generators implementation +////////////////////////////////////////////////////////////////////// +/* #define TRACERPLUS (-5) */ +#include "lasp_siggen_impl.h" +#include "debugtrace.hpp" +#include "lasp_mathtypes.h" + +DEBUGTRACE_VARIABLES; + +/** The fixed number of Newton iterations t.b.d. for tuning the sweep start and + * stop frequency in logarithmic sweeps */ +#define NITER_NEWTON 20 + +Noise::Noise(){DEBUGTRACE_ENTER} + +vd Noise::genSignalUnscaled(us nframes) { + return arma::randn(nframes); +} +void Noise::resetImpl() {} + +Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) {} + +vd Sine::genSignalUnscaled(const us nframes) { + const d pi = arma::datum::pi; + vd phase_vec = + arma::linspace(phase, phase + omg * (nframes - 1) / fs, nframes); + phase += omg * nframes / fs; + while (phase > 2 * arma::datum::pi) { + phase -= 2 * pi; + } + return arma::sin(phase_vec); +} + +vd Periodic::genSignalUnscaled(const us nframes) { + us signal_idx = 0; + vd res(nframes); + while (signal_idx < nframes) { + res[signal_idx] = _signal[_cur_pos]; + _cur_pos++; + _cur_pos %= _signal.size(); + } + return res; +} + +Sweep::Sweep(const d fl, const d fu, const d Ts, const d Tq, const us flags) + : fl_(fl), fu_(fu), Ts(Ts), Tq(Tq), flags(flags) { + if (fl <= 0 || fu < fl || Ts <= 0) { + throw std::runtime_error("Invalid sweep parameters"); + } + if ((flags & ForwardSweep) && (flags & BackwardSweep)) { + throw std::runtime_error( + "Both forward and backward sweep flag set. Please only set either one " + "or none for a continuous sweep"); + } +} + +void Sweep::resetImpl() { + + _cur_pos = 0; + + bool forward_sweep = flags & ForwardSweep; + bool backward_sweep = flags & BackwardSweep; + + const d Dt = 1 / fs; // Deltat + + // Estimate N, the number of samples in the sweep part (non-quiescent part): + const us Ns = (us)(Ts * fs); + const us Nq = (us)(Tq * fs); + const us N = Ns + Nq; + + _signal = vd(N, arma::fill::zeros); + index = 0; + + d fl, fu; + /* Swap fl and fu for a backward sweep */ + if (backward_sweep) { + fu = fl_; + fl = fu_; + } else { + /* Case of continuous sweep, or forward sweep */ + fl = fl_; + fu = fu_; + } + + d phase = 0; + + /* Linear sweep */ + if (flags & LinearSweep) { + if (forward_sweep || backward_sweep) { + /* Forward or backward sweep */ + /* TRACE(15, "Forward or backward sweep"); */ + us K = (us)(Dt * (fl * N + 0.5 * (N - 1) * (fu - fl))); + d eps_num = ((d)K) / Dt - fl * N - 0.5 * (N - 1) * (fu - fl); + d eps = eps_num / (0.5 * (N - 1)); + /* iVARTRACE(15, K); */ + /* dVARTRACE(15, eps); */ + + for (us n = 0; n < Ns; n++) { + _signal(n) = d_sin(phase); + d fn = fl + ((d)n) / N * (fu + eps - fl); + phase += 2 * arma::datum::pi * Dt * fn; + } + } else { + /* Continous sweep */ + /* TRACE(15, "continuous sweep"); */ + + /* iVARTRACE(17, N); */ + /* dVARTRACE(17, fl); */ + /* dVARTRACE(17, fu); */ + + const us Nf = Ns / 2; + const us Nb = Ns - Nf; + + /* Phi halfway */ + d phih = 2 * number_pi * Dt * (fl * Nf + 0.5 * (Nf - 1) * (fu - fl)); + + us K = + (us)(phih / (2 * number_pi) + Dt * (fu * Nb - (Nb - 1) * (fu - fl))); + + d eps_num1 = (K - phih / (2 * number_pi)) / Dt; + d eps_num2 = -fu * Nb + (Nb - 1) * (fu - fl); + d eps = (eps_num1 + eps_num2) / (0.5 * (Nb + 1)); + + /* iVARTRACE(15, K); */ + /* dVARTRACE(15, eps); */ + d phase = 0; + + for (us n = 0; n <= Ns; n++) { + /* iVARTRACE(17, n); */ + if (n < N) { + _signal[n] = d_sin(phase); + } + + d fn; + if (n <= Nf) { + fn = fl + ((d)n) / Nf * (fu - fl); + } else { + fn = fu - ((d)n - Nf) / Nb * (fu + eps - fl); + } + /* dbgassert(fn >= 0, "BUG"); */ + + phase += 2 * number_pi * Dt * fn; + } + /* This should be a very small number!! */ + /* dVARTRACE(15, phase); */ + } + } else if (flags & LogSweep) { + + DEBUGTRACE_PRINT("Exponential sweep"); + if (forward_sweep || backward_sweep) { + /* Forward or backward sweep */ + DEBUGTRACE_PRINT("Forward or backward sweep"); + d k1 = (fu / fl); + us K = (us)(Dt * fl * (k1 - 1) / (d_pow(k1, 1.0 / N) - 1)); + d k = k1; + + /* Iterate k to the right solution */ + d E; + for (us iter = 0; iter < 10; iter++) { + E = 1 + K / (Dt * fl) * (d_pow(k, 1.0 / N) - 1) - k; + d dEdk = K / (Dt * fl) * d_pow(k, 1.0 / N) / (N * k) - 1; + k -= E / dEdk; + } + + DEBUGTRACE_PRINT(K); + DEBUGTRACE_PRINT(K1); + DEBUGTRACE_PRINT(k); + DEBUGTRACE_PRINT(E); + + for (us n = 0; n < Ns; n++) { + _signal[n] = d_sin(phase); + d fn = fl * d_pow(k, ((d)n) / N); + phase += 2 * number_pi * Dt * fn; + } + } else { + + DEBUGTRACE_PRINT("Continuous sweep"); + + const us Nf = N / 2; + const us Nb = N - Nf; + const d k1 = (fu / fl); + const d phif1 = + 2 * number_pi * Dt * fl * (k1 - 1) / (d_pow(k1, 1.0 / Nf) - 1); + const us K = (us)(phif1 / (2 * number_pi) + + Dt * fu * (1 / k1 - 1) / (d_pow(1 / k1, 1.0 / Nb) - 1)); + + d E; + d k = k1; + + /* Newton iterations to converge k to the value such that the sweep is + * continuous */ + for (us iter = 0; iter < NITER_NEWTON; iter++) { + E = (k - 1) / (d_pow(k, 1.0 / Nf) - 1) + + (k - 1) / (1 - d_pow(k, -1.0 / Nb)) - K / Dt / fl; + DEBUGTRACE_PRINT(E); + + /* All parts of the derivative of above error E to k */ + d dEdk1 = 1 / (d_pow(k, 1.0 / Nf) - 1); + d dEdk2 = (1 / k - 1) / (d_pow(k, -1.0 / Nb) - 1); + d dEdk3 = -1 / (k * (d_pow(k, -1.0 / Nb) - 1)); + d dEdk4 = d_pow(k, -1.0 / Nb) * (1 / k - 1) / + (Nb * d_pow(d_pow(k, -1.0 / Nb) - 1, 2)); + d dEdk5 = -d_pow(k, 1.0 / Nf) * (k - 1) / + (Nf * k * d_pow(d_pow(k, 1.0 / Nf) - 1, 2)); + d dEdk = dEdk1 + dEdk2 + dEdk3 + dEdk4 + dEdk5; + + /* Iterate! */ + k -= E / dEdk; + + } + + DEBUGTRACE_PRINT(K); + DEBUGTRACE_PRINT(K1); + DEBUGTRACE_PRINT(k); + DEBUGTRACE_PRINT(E); + + for (us n = 0; n <= Ns; n++) { + /* iVARTRACE(17, n); */ + if (n < Ns) { + _signal[n] = d_sin(phase); + } + + d fn; + if (n <= Nf) { + fn = fl * d_pow(k, ((d)n) / Nf); + } else { + fn = fl * k * d_pow(1 / k, ((d)n - Nf) / Nb); + } + /* dbgassert(fn >= 0, "BUG"); */ + + phase += 2 * number_pi * Dt * fn; + while (phase > 2 * number_pi) + phase -= 2 * number_pi; + /* dVARTRACE(17, phase); */ + } + /* This should be a very small number!! */ + DEBUGTRACE_PRINT(phase); + } + } +} diff --git a/lasp/dsp/lasp_siggen_impl.h b/lasp/dsp/lasp_siggen_impl.h new file mode 100644 index 0000000..9bc38ff --- /dev/null +++ b/lasp/dsp/lasp_siggen_impl.h @@ -0,0 +1,79 @@ +#pragma once +#include "lasp_siggen.h" +#include "lasp_types.h" + +class Noise : public Siggen { + d level_linear; + virtual vd genSignalUnscaled(const us nframes) override; + void resetImpl() override; + public: + + /** + * @brief Constructs a noise generator. If no filter is used, the output will + * be white noise. By default, the output will be standard deviation = 1 + * noise, which clips the output for standard audio devices, so make sure the + * level is set properly. + */ + Noise(); + ~Noise() = default; + +}; + +class Sine : public Siggen { + d phase = 0; + d omg; + + public: + + /** + * @brief Create a sine wave generator + * + * @param freq_Hz + * @param level_dB + */ + Sine(const d freq_Hz); + ~Sine() = default; + virtual vd genSignalUnscaled(const us nframes) override; + void setFreq(const d newFreq); +}; +class Periodic: public Siggen { + protected: + vd _signal { 1, arma::fill::zeros}; + us _cur_pos = 0; + public: + + virtual vd genSignalUnscaled(const us nframes) override; + ~Periodic() = default; + +}; + +class Sweep : public Periodic { + d fl_, fu_, Ts, Tq; + us index; + us flags; + + static constexpr int ForwardSweep = 1 << 0; + static constexpr int BackwardSweep = 1 << 1; + static constexpr int LinearSweep = 1 << 2; + static constexpr int LogSweep = 1 << 3; + + void resetImpl() override; + + public: + + /** + * Create a sweep signal + * + * @param[in] fl: Lower frequency [Hz] + * @param[in] fu: Upper frequency [Hz] + * @param[in] Ts: Sweep time [s] + * @param[in] Tq: Quescent tail time [s]. Choose this value long enough to + * avoid temporal aliasing in case of measuring impulse responses. + * @param[in] sweep_flags: Sweep period [s] + */ + Sweep(const d fl, const d fu, const d Ts, const d Tq, + const us sweep_flags); + + ~Sweep() = default; + +}; diff --git a/lasp/c/lasp_types.h b/lasp/dsp/lasp_types.h similarity index 87% rename from lasp/c/lasp_types.h rename to lasp/dsp/lasp_types.h index b937eb3..c350878 100644 --- a/lasp/c/lasp_types.h +++ b/lasp/dsp/lasp_types.h @@ -51,14 +51,6 @@ typedef double complex c; #endif -/// I need these numbers so often, that they can be in the global -/// namespace. -#define LASP_SUCCESS 0 -#define LASP_INTERRUPTED (-3) -#define LASP_MALLOC_FAILED (-1) -#define LASP_FAILURE (-2) - - #endif // LASP_TYPES_H ////////////////////////////////////////////////////////////////////// diff --git a/lasp/dsp/lasp_window.cpp b/lasp/dsp/lasp_window.cpp new file mode 100644 index 0000000..50a5860 --- /dev/null +++ b/lasp/dsp/lasp_window.cpp @@ -0,0 +1,57 @@ +// lasp_window.cpp +// +// Author: J.A. de Jong - ASCEE +// +#include "lasp_window.h" + +// Safe some typing. Linspace form 0 up to (and NOT including N). +#define lin0N arma::linspace(0, N - 1, N) + +vd Window::hann(const us N) { + return arma::pow(arma::sin(number_pi * lin0N), 2); +} +vd Window::hamming(const us N) { + d alpha = 25.0 / 46.0; + return alpha - (1 - alpha) * arma::cos(2 * number_pi * lin0N / (N - 1)); +} +vd Window::blackman(const us N) { + d a0 = 7938. / 18608.; + d a1 = 9240. / 18608.; + d a2 = 1430. / 18608.; + return a0 - a1 * d_cos(2 * number_pi * lin0N / (N - 1)) + + a2 * d_cos(4 * number_pi * lin0N / (N - 1)); +} + +vd Window::rectangular(const us N) { return arma::ones(N); } + +vd Window::bartlett(const us N) { + return 1 - arma::abs(2 * (lin0N - (N - 1) / 2.) / (N - 1)); +} +vd Window::create(const WindowType w, const us N) { + + switch (w) { + case WindowType::Hann: { + return hann(N); + break; + } + case WindowType::Hamming: { + return hamming(N); + break; + } + case WindowType::Rectangular: { + return rectangular(N); + break; + } + case WindowType::Bartlett: { + return bartlett(N); + break; + } + case WindowType::Blackman: { + return blackman(N); + break; + } + default: + abort(); + break; + } +} diff --git a/lasp/dsp/lasp_window.h b/lasp/dsp/lasp_window.h new file mode 100644 index 0000000..514dbfe --- /dev/null +++ b/lasp/dsp/lasp_window.h @@ -0,0 +1,73 @@ +#pragma once +#include "lasp_mathtypes.h" + +/** + * @brief Window (aka taper) functions of a certain type + */ +class Window { + + public: + enum class WindowType { + Hann = 0, + Hamming = 1, + Rectangular = 2, + Bartlett = 3, + Blackman = 4, + + }; + + /** + * @brief Dispatcher: create a window based on enum type and len + * + * @param w Window type + * @param len Length of the window (typically, integer power of two). + * + * @return Window vector of values + */ + static vd create(const WindowType w,const us len); + + /** + * @brief Hann window + * + * @param len Length of the window (typically, integer power of two). + * + * @return vector of values + */ + static vd hann(const us len); + + /** + * @brief Hamming window + * + * @param len Length of the window (typically, integer power of two). + * + * @return vector of values + */ + static vd hamming(const us len); + /** + * @brief Rectangular (boxcar) window. + * + * @param len Length of the window (typically, integer power of two). + * + * @return vector of values + */ + static vd rectangular(const us len); + + /** + * @brief Bartlett window. + * + * @param len Length of the window (typically, integer power of two). + * + * @return vector of values + */ + static vd bartlett(const us len); + + /** + * @brief Blackman window. + * + * @param len Length of the window (typically, integer power of two). + * + * @return vector of values + */ + static vd blackman(const us len); + +}; diff --git a/lasp/lasp_common.py b/lasp/lasp_common.py index 3ce3c3d..38c300a 100644 --- a/lasp/lasp_common.py +++ b/lasp/lasp_common.py @@ -2,11 +2,11 @@ import shelve, logging, sys, appdirs, os, platform import numpy as np -from .wrappers import Window as wWindow from collections import namedtuple from dataclasses import dataclass from dataclasses_json import dataclass_json from enum import Enum, unique, auto +from .dsp import Window as wWindow """ Common definitions used throughout the code. @@ -271,11 +271,11 @@ class this_lasp_shelve(Shelve): @unique class Window(Enum): - hann = (wWindow.hann, 'Hann') - hamming = (wWindow.hamming, 'Hamming') - rectangular = (wWindow.rectangular, 'Rectangular') - bartlett = (wWindow.bartlett, 'Bartlett') - blackman = (wWindow.blackman, 'Blackman') + hann = (wWindow.Hann, 'Hann') + hamming = (wWindow.Hamming, 'Hamming') + rectangular = (wWindow.Rectangular, 'Rectangular') + bartlett = (wWindow.Bartlett, 'Bartlett') + blackman = (wWindow.Blackman, 'Blackman') @staticmethod def fillComboBox(cb): diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx deleted file mode 100644 index 444bf49..0000000 --- a/lasp/wrappers.pyx +++ /dev/null @@ -1,798 +0,0 @@ -""" -This file contains the Cython wrapper functions to C implementations. -""" -include "config.pxi" -from libc.stdio cimport printf -from numpy cimport import_array - -import_array() - - -IF LASP_DOUBLE_PRECISION == "ON": - ctypedef double d - ctypedef double complex c - NUMPY_FLOAT_TYPE = np.float64 - NUMPY_COMPLEX_TYPE = np.complex128 - CYTHON_NUMPY_FLOAT_t = cnp.NPY_FLOAT64 - CYTHON_NUMPY_COMPLEX_t = cnp.NPY_COMPLEX128 - -ELSE: - ctypedef float d - ctypedef float complex c - NUMPY_FLOAT_TYPE = np.float32 - NUMPY_COMPLEX_TYPE = np.complex64 - CYTHON_NUMPY_FLOAT_t = cnp.NPY_FLOAT32 - CYTHON_NUMPY_COMPLEX_t = cnp.NPY_COMPLEX64 - -ctypedef size_t us - -cdef extern from "lasp_tracer.h": - void setTracerLevel(int) - void TRACE(int,const char*) - void fsTRACE(int) - void feTRACE(int) - void clearScreen() - - -cdef extern from "lasp_mat.h" nogil: - ctypedef struct dmat: - us n_cols - us n_rows - d* _data - bint _foreign_data - ctypedef struct cmat: - pass - ctypedef cmat vc - ctypedef dmat vd - - dmat dmat_foreign_data(us n_rows, - us n_cols, - d* data, - bint own_data) - cmat cmat_foreign_data(us n_rows, - us n_cols, - c* data, - bint own_data) - - cmat cmat_alloc(us n_rows,us n_cols) - dmat dmat_alloc(us n_rows,us n_cols) - vd vd_foreign(const us size,d* data) - void vd_free(vd*) - - void dmat_free(dmat*) - void cmat_free(cmat*) - void cmat_copy(cmat* to,cmat* from_) - - -cdef extern from "numpy/arrayobject.h": - void PyArray_ENABLEFLAGS(cnp.ndarray arr, int flags) - - -cdef extern from "lasp_python.h": - object dmat_to_ndarray(dmat*,bint transfer_ownership) - -__all__ = ['AvPowerSpectra', 'SosFilterBank', 'FilterBank', 'Siggen', - 'sweep_flag_forward', 'sweep_flag_backward', 'sweep_flag_linear', - 'sweep_flag_exponential', - 'load_fft_wisdom', 'store_fft_wisdom'] - -setTracerLevel(15) - -cdef extern from "cblas.h": - int openblas_get_num_threads() - void openblas_set_num_threads(int) - -# If we touch this variable: we get segfaults when running from -# Spyder! -# openblas_set_num_threads(8) -# print("Number of threads: ", -# openblas_get_num_threads()) - -def cls(): - clearScreen() -# cls() - -cdef extern from "lasp_fft.h": - void c_load_fft_wisdom "load_fft_wisdom" (const char* wisdom) - char* c_store_fft_wisdom "store_fft_wisdom" () - ctypedef struct c_Fft "Fft" - c_Fft* Fft_create(us nfft) - void Fft_free(c_Fft*) - void Fft_fft(c_Fft*,dmat * timedate,cmat * res) nogil - void Fft_ifft(c_Fft*,cmat * freqdata,dmat* timedata) nogil - us Fft_nfft(c_Fft*) - -def load_fft_wisdom(const unsigned char[::1] wisdom): - c_load_fft_wisdom( &wisdom[0]) - -from cpython cimport PyBytes_FromString -from libc.stdlib cimport free - -def store_fft_wisdom(): - cdef char* wisdom = c_store_fft_wisdom() - - if wisdom != NULL: - try: - bts = PyBytes_FromString(wisdom) - finally: - free(wisdom) - return bts - else: - return None - -cdef class Fft: - cdef: - c_Fft* _fft - - def __cinit__(self, us nfft): - self._fft = Fft_create(nfft) - if self._fft == NULL: - raise RuntimeError('Fft allocation failed') - - def __dealloc__(self): - if self._fft!=NULL: - Fft_free(self._fft) - - def fft(self,d[::1,:] timedata): - - cdef us nfft = Fft_nfft(self._fft) - cdef us nchannels = timedata.shape[1] - assert timedata.shape[0] ==nfft - - result = np.empty((nfft//2+1,nchannels), - dtype=NUMPY_COMPLEX_TYPE, - order='F') - - # result[:,:] = np.nan+1j*np.nan - cdef c[::1,:] result_view = result - cdef cmat r = cmat_foreign_data(result.shape[0], - result.shape[1], - &result_view[0,0], - False) - - cdef dmat t = dmat_foreign_data(timedata.shape[0], - timedata.shape[1], - &timedata[0,0], - False) - with nogil: - Fft_fft(self._fft,&t,&r) - - dmat_free(&t) - cmat_free(&r) - - return result - - def ifft(self,c[::1,:] freqdata): - - cdef us nfft = Fft_nfft(self._fft) - cdef us nchannels = freqdata.shape[1] - assert freqdata.shape[0] == nfft//2+1 - - - # result[:,:] = np.nan+1j*np.nan - - cdef cmat f = cmat_foreign_data(freqdata.shape[0], - freqdata.shape[1], - &freqdata[0,0], - False) - - timedata = np.empty((nfft,nchannels), - dtype=NUMPY_FLOAT_TYPE, - order='F') - - cdef d[::1,:] timedata_view = timedata - cdef dmat t = dmat_foreign_data(timedata.shape[0], - timedata.shape[1], - &timedata_view[0,0], - False) - - with nogil: - Fft_ifft(self._fft,&f,&t) - - dmat_free(&t) - cmat_free(&f) - - return timedata - -cdef extern from "lasp_window.h": - ctypedef enum WindowType: - Hann - Hamming - Rectangular - Bartlett - Blackman - -# Export these constants to Python -class Window: - hann = Hann - hamming = Hamming - rectangular = Rectangular - bartlett = Bartlett - blackman = Blackman - -cdef extern from "lasp_ps.h": - ctypedef struct c_PowerSpectra "PowerSpectra" - c_PowerSpectra* PowerSpectra_alloc(const us nfft, - const WindowType wt) - - void PowerSpectra_compute(const c_PowerSpectra* ps, - const dmat * timedata, - cmat * result) nogil - - - void PowerSpectra_free(c_PowerSpectra*) - -cdef class PowerSpectra: - cdef: - c_PowerSpectra* _ps - - def __cinit__(self, us nfft,us window=Window.rectangular): - self._ps = PowerSpectra_alloc(nfft, window) - if self._ps == NULL: - raise RuntimeError('PowerSpectra allocation failed') - - def compute(self,d[::1,:] timedata): - cdef: - us nchannels = timedata.shape[1] - us nfft = timedata.shape[0] - int rv - dmat td - cmat result_mat - - - td = dmat_foreign_data(nfft, - nchannels, - &timedata[0,0], - False) - - # The array here is created in such a way that the strides - # increase with increasing dimension. This is required for - # interoperability with the C-code, that stores all - # cross-spectra in a 2D matrix, where the first axis is the - # frequency axis, and the second axis corresponds to a certain - # cross-spectrum, as C_ij(f) = result[freq,i+j*nchannels] - - result = np.empty((nfft//2+1,nchannels,nchannels), - dtype = NUMPY_COMPLEX_TYPE, - order='F') - - cdef c[::1,:,:] result_view = result - - result_mat = cmat_foreign_data(nfft//2+1, - nchannels*nchannels, - &result_view[0,0,0], - False) - - - with nogil: - PowerSpectra_compute(self._ps,&td,&result_mat) - - dmat_free(&td) - cmat_free(&result_mat) - - return result - - - def __dealloc__(self): - if self._ps != NULL: - PowerSpectra_free(self._ps) - - -cdef extern from "lasp_aps.h": - ctypedef struct c_AvPowerSpectra "AvPowerSpectra" - c_AvPowerSpectra* AvPowerSpectra_alloc(const us nfft, - const us nchannels, - d overlap_percentage, - const WindowType wt, - const vd* weighting) - - cmat* AvPowerSpectra_addTimeData(const c_AvPowerSpectra* ps, - const dmat * timedata) nogil - - - void AvPowerSpectra_free(c_AvPowerSpectra*) - us AvPowerSpectra_getAverages(const c_AvPowerSpectra*); - -cdef class AvPowerSpectra: - cdef: - c_AvPowerSpectra* aps - us nfft, nchannels - - def __cinit__(self,us nfft, - us nchannels, - d overlap_percentage, - us window=Window.hann, - d[:] weighting = np.array([])): - - - cdef vd weighting_vd - cdef vd* weighting_ptr = NULL - if(weighting.size != 0): - weighting_vd = dmat_foreign_data(weighting.size,1, - &weighting[0],False) - weighting_ptr = &weighting_vd - - self.aps = AvPowerSpectra_alloc(nfft, - nchannels, - overlap_percentage, - window, - weighting_ptr) - self.nchannels = nchannels - self.nfft = nfft - - if self.aps == NULL: - raise RuntimeError('AvPowerSpectra allocation failed') - - def __dealloc__(self): - if self.aps: - AvPowerSpectra_free(self.aps) - def getAverages(self): - return AvPowerSpectra_getAverages(self.aps) - - def addTimeData(self,d[::1,:] timedata): - """! - Adds time data, returns current result - """ - cdef: - us nsamples = timedata.shape[0] - us nchannels = timedata.shape[1] - dmat td - cmat* result_ptr - - if nchannels != self.nchannels: - raise RuntimeError('Invalid number of channels') - - td = dmat_foreign_data(nsamples, - nchannels, - &timedata[0,0], - False) - - result = np.empty((self.nfft//2+1,nchannels,nchannels), - dtype = NUMPY_COMPLEX_TYPE, - order='F') - - cdef c[::1,:,:] result_view = result - - cdef cmat res = cmat_foreign_data(self.nfft//2+1, - nchannels*nchannels, - &result_view[0,0,0], - False) - with nogil: - result_ptr = AvPowerSpectra_addTimeData(self.aps, - &td) - - # The array here is created in such a way that the strides - # increase with increasing dimension. This is required for - # interoperability with the C-code, that stores all - # cross-spectra in a 2D matrix, where the first axis is the - # frequency axis, and the second axis corresponds to a certain - # cross-spectrum, as C_ij(f) = result[freq,i+j*nchannels] - - # Copy result - cmat_copy(&res,result_ptr) - - cmat_free(&res) - dmat_free(&td) - - return result - -cdef extern from "lasp_firfilterbank.h": - ctypedef struct c_Firfilterbank "Firfilterbank" - c_Firfilterbank* Firfilterbank_create(const dmat* h,const us nfft) nogil - dmat Firfilterbank_filter(c_Firfilterbank* fb,const vd* x) nogil - void Firfilterbank_free(c_Firfilterbank* fb) nogil - - -cdef class FilterBank: - cdef: - c_Firfilterbank* fb - def __cinit__(self,d[::1,:] h, us nfft): - cdef dmat hmat = dmat_foreign_data(h.shape[0], - h.shape[1], - &h[0,0], - False) - - self.fb = Firfilterbank_create(&hmat,nfft) - dmat_free(&hmat) - if not self.fb: - raise RuntimeError('Error creating FilberBank') - - def __dealloc__(self): - if self.fb: - Firfilterbank_free(self.fb) - - def filter_(self,d[::1, :] input_): - assert input_.shape[1] == 1 - cdef dmat input_vd = dmat_foreign_data(input_.shape[0],1, - &input_[0, 0],False) - - - cdef dmat output - with nogil: - output = Firfilterbank_filter(self.fb,&input_vd) - - # Steal the pointer from output - result = dmat_to_ndarray(&output,True) - - dmat_free(&output) - vd_free(&input_vd) - - return result - -cdef extern from "lasp_sosfilterbank.h": - ctypedef struct c_Sosfilterbank "Sosfilterbank" - c_Sosfilterbank* Sosfilterbank_create(const us nthreads, const us filterbank_size, const us nsections) nogil - void Sosfilterbank_setFilter(c_Sosfilterbank* fb,const us filter_no, const vd filter_coefs) - dmat Sosfilterbank_filter(c_Sosfilterbank* fb,const vd* x) nogil - void Sosfilterbank_free(c_Sosfilterbank* fb) nogil - - -cdef class SosFilterBank: - cdef: - c_Sosfilterbank* fb - us nsections - def __cinit__(self,const us filterbank_size, const us nsections): - self.nsections = nsections - - self.fb = Sosfilterbank_create(0, filterbank_size,nsections) - - def setFilter(self,us filter_no, d[:, ::1] sos): - """ - - Args: - filter_no: Filter number of the filterbank to set the - filter for - sos: Second axis are the filter coefficients, first axis - is the section, second axis are the coefficients for that section. - Storage is in agreement with specification from Scipy: first axis - is the section, second axis are the coefficients for that section. - - """ - if sos.shape[0] != self.nsections: - raise RuntimeError(f'Invalid number of sections in filter data, should be {self.nsections}.') - elif sos.shape[1] != 6: - raise RuntimeError('Illegal number of filter coefficients in section. Should be 6.') - cdef dmat coefs = dmat_foreign_data(sos.size,1, - &sos[0, 0],False # No copying - ) - Sosfilterbank_setFilter(self.fb,filter_no, coefs) - - def __dealloc__(self): - if self.fb: - Sosfilterbank_free(self.fb) - - def filter_(self,d[::1, :] input_): - # Only single channel input - assert input_.shape[1] == 1 - cdef dmat input_vd = dmat_foreign_data(input_.shape[0],1, - &input_[0, 0],False) - - - cdef dmat output - with nogil: - output = Sosfilterbank_filter(self.fb,&input_vd) - #printf('Came back from filter\n') - # Steal the pointer from output - result = dmat_to_ndarray(&output,True) - #printf('Converted to array\n') - dmat_free(&output) - vd_free(&input_vd) - #printf('Ready to return\n') - return result - -cdef extern from "lasp_decimation.h": - ctypedef struct c_Decimator "Decimator" - ctypedef enum DEC_FAC: - DEC_FAC_4 - - c_Decimator* Decimator_create(us nchannels,DEC_FAC d) nogil - dmat Decimator_decimate(c_Decimator* dec,const dmat* samples) nogil - void Decimator_free(c_Decimator* dec) nogil - - -cdef extern from "lasp_slm.h" nogil: - ctypedef struct c_Slm "Slm" - d TAU_FAST, TAU_SLOW, TAU_IMPULSE - c_Slm* Slm_create(c_Sosfilterbank* prefilter, - c_Sosfilterbank* bandpass, - d fs, d tau, d ref_level, - us* downsampling_fac) - dmat Slm_run(c_Slm* slm,vd* input_data) - void Slm_free(c_Slm* slm) - vd Slm_Leq(c_Slm*) - vd Slm_Lmax(c_Slm*) - vd Slm_Lpeak(c_Slm*) - - -tau_fast = TAU_FAST -tau_slow = TAU_SLOW -tau_impulse = TAU_IMPULSE - -cdef class Slm: - cdef: - c_Slm* c_slm - public us downsampling_fac - - def __cinit__(self, d[::1] sos_prefilter, - d[:, ::1] sos_bandpass, - d fs, d tau, d ref_level): - - cdef: - us prefilter_nsections - us bandpass_nsections - us bandpass_nchannels - c_Sosfilterbank* prefilter = NULL - c_Sosfilterbank* bandpass = NULL - vd coefs_vd - d[:] coefs - - if sos_prefilter is not None: - assert sos_prefilter.size % 6 == 0 - prefilter_nsections = sos_prefilter.size // 6 - prefilter = Sosfilterbank_create(0, 1,prefilter_nsections) - coefs = sos_prefilter - coefs_vd = dmat_foreign_data(prefilter_nsections*6,1, - &coefs[0],False) - Sosfilterbank_setFilter(prefilter, 0, coefs_vd) - - if prefilter is NULL: - raise RuntimeError('Error creating pre-filter') - - if sos_bandpass is not None: - assert sos_bandpass.shape[1] % 6 == 0 - bandpass_nsections = sos_bandpass.shape[1] // 6 - bandpass_nchannels = sos_bandpass.shape[0] - bandpass = Sosfilterbank_create(0, - bandpass_nchannels, - bandpass_nsections) - if bandpass == NULL: - if prefilter: - Sosfilterbank_free(prefilter) - raise RuntimeError('Error creating bandpass filter') - - - for i in range(bandpass_nchannels): - coefs = sos_bandpass[i, :] - coefs_vd = dmat_foreign_data(bandpass_nsections*6,1, - &coefs[0],False) - - Sosfilterbank_setFilter(bandpass, i, coefs_vd) - - self.c_slm = Slm_create(prefilter, bandpass, - fs, tau, ref_level, - &self.downsampling_fac) - if self.c_slm is NULL: - Sosfilterbank_free(prefilter) - Sosfilterbank_free(bandpass) - raise RuntimeError('Error creating sound level meter') - - def run(self, d[:, ::1] data): - assert data.shape[1] == 1 - cdef vd data_vd = dmat_foreign_data(data.shape[0], 1, - &data[0,0], False) - cdef dmat res - with nogil: - res = Slm_run(self.c_slm, &data_vd) - result = dmat_to_ndarray(&res,True) - return result - - def Leq(self): - cdef vd res - res = Slm_Leq(self.c_slm) - # True below, means transfer ownership of allocated data to Numpy - return dmat_to_ndarray(&res,True) - - def Lmax(self): - cdef vd res - res = Slm_Lmax(self.c_slm) - # True below, means transfer ownership of allocated data to Numpy - return dmat_to_ndarray(&res,True) - - def Lpeak(self): - cdef vd res - res = Slm_Lpeak(self.c_slm) - # True below, means transfer ownership of allocated data to Numpy - return dmat_to_ndarray(&res,True) - - def __dealloc__(self): - if self.c_slm: - Slm_free(self.c_slm) - - - - -cdef class Decimator: - cdef: - c_Decimator* dec - us nchannels - def __cinit__(self, us nchannels,us dec_fac): - assert dec_fac == 4, 'Invalid decimation factor' - self.nchannels = nchannels - self.dec = Decimator_create(nchannels,DEC_FAC_4) - if not self.dec: - raise RuntimeError('Error creating decimator') - - def decimate(self,d[::1,:] samples): - assert samples.shape[1] == self.nchannels,'Invalid number of channels' - if samples.shape[0] == 0: - return np.zeros((0, self.nchannels)) - - cdef dmat d_samples = dmat_foreign_data(samples.shape[0], - samples.shape[1], - &samples[0,0], - False) - - cdef dmat res = Decimator_decimate(self.dec,&d_samples) - result = dmat_to_ndarray(&res,True) - dmat_free(&res) - return result - - def __dealloc__(self): - if self.dec != NULL: - Decimator_free(self.dec) - - -cdef extern from "lasp_siggen.h": - ctypedef struct c_Siggen "Siggen" - us SWEEP_FLAG_FORWARD, SWEEP_FLAG_BACKWARD, SWEEP_FLAG_LINEAR - us SWEEP_FLAG_EXPONENTIAL - - c_Siggen* Siggen_Noise_create(d fs, d level_dB, c_Sosfilterbank* - colorfilter) - c_Siggen* Siggen_Sinewave_create(d fs, d freq, d level_dB) - c_Siggen* Siggen_Sweep_create(d fs, d fl, - d fu, d Ts,d Tq, us sweep_flags, - d level_dB) - void Siggen_setLevel(c_Siggen*,d new_level_dB) - us Siggen_getN(const c_Siggen*) - void Siggen_genSignal(c_Siggen*, vd* samples) nogil - void Siggen_free(c_Siggen*) - -# Sweep flags -sweep_flag_forward = SWEEP_FLAG_FORWARD -sweep_flag_backward = SWEEP_FLAG_BACKWARD - -sweep_flag_linear = SWEEP_FLAG_LINEAR -sweep_flag_exponential = SWEEP_FLAG_EXPONENTIAL - -from .filter import PinkNoise - -cdef class Siggen: - - cdef c_Siggen *_siggen - - def __cinit__(self): - self._siggen = NULL - - def __dealloc__(self): - if self._siggen: - Siggen_free(self._siggen) - - def setLevel(self,d level_dB): - Siggen_setLevel(self._siggen, level_dB) - - def genSignal(self, us nsamples): - output = np.empty(nsamples, dtype=np.float) - assert self._siggen != NULL - - cdef d[:] output_view = output - - cdef dmat output_dmat = dmat_foreign_data(nsamples, - 1, - &output_view[0], - False) - with nogil: - Siggen_genSignal(self._siggen, - &output_dmat) - - return output - - def getN(self): - return Siggen_getN(self._siggen) - - def progress(self): - """ - TODO: Should be implemented to return the current position in the - generator. - - """ - return None - - @staticmethod - def sineWave(d fs,d freq,d level_dB): - cdef c_Siggen* c_siggen = Siggen_Sinewave_create(fs, freq, level_dB) - siggen = Siggen() - siggen._siggen = c_siggen - return siggen - - - @staticmethod - def noise(d fs, d level_dB, d[::1] colorfilter_coefs=None): - cdef: - c_Sosfilterbank* colorfilter = NULL - if colorfilter_coefs is not None: - assert colorfilter_coefs.size % 6 == 0 - colorfilter_nsections = colorfilter_coefs.size // 6 - colorfilter = Sosfilterbank_create(0, 1,colorfilter_nsections) - coefs = colorfilter_coefs - coefs_vd = dmat_foreign_data(colorfilter_nsections*6,1, - &colorfilter_coefs[0],False) - Sosfilterbank_setFilter(colorfilter, 0, coefs_vd) - - if colorfilter is NULL: - raise RuntimeError('Error creating pre-filter') - cdef c_Siggen* c_siggen = Siggen_Noise_create(fs, level_dB, colorfilter) - siggen = Siggen() - siggen._siggen = c_siggen - return siggen - - - @staticmethod - def sweep(d fs, d fl, d fu, d Ts, d Tq, us sweep_flags, d level_dB): - cdef c_Siggen* c_siggen = Siggen_Sweep_create(fs, - fl, - fu, - Ts, - Tq, - sweep_flags, - level_dB) - if c_siggen == NULL: - raise ValueError('Failed creating signal generator') - siggen = Siggen() - siggen._siggen = c_siggen - return siggen - -cdef extern from "lasp_eq.h" nogil: - ctypedef struct c_Eq "Eq" - c_Eq* Eq_create(c_Sosfilterbank* fb) - vd Eq_equalize(c_Eq* eq,const vd* input_data) - us Eq_getNLevels(const c_Eq* eq) - void Eq_setLevels(c_Eq* eq, const vd* levels) - void Eq_free(c_Eq* eq) - -cdef class Equalizer: - cdef: - c_Eq* ceq - - def __cinit__(self, SosFilterBank cdef_fb): - """ - Initialize equalizer using given filterbank. Note: Steals pointer of - underlying c_Sosfilterbank!! - """ - self.ceq = Eq_create(cdef_fb.fb) - # Set this pointer to NULL, such that the underlying c_SosfilterBank is - # not deallocated. - cdef_fb.fb = NULL - - def getNLevels(self): - return Eq_getNLevels(self.ceq) - - def setLevels(self,d[:] new_levels): - cdef dmat dmat_new_levels = dmat_foreign_data(new_levels.shape[0], - 1, - &new_levels[0], - False) - Eq_setLevels(self.ceq, &dmat_new_levels) - dmat_free(&dmat_new_levels) - - def equalize(self, d[::1] input_data): - cdef: - vd res - cdef dmat input_dmat = dmat_foreign_data(input_data.size, - 1, - &input_data[0], - False) - with nogil: - res = Eq_equalize(self.ceq, &input_dmat) - - # Steal the pointer from output - py_res = dmat_to_ndarray(&res,True)[:,0] - - dmat_free(&res) - vd_free(&input_dmat) - - return py_res - - def __dealloc__(self): - if self.ceq: - Eq_free(self.ceq) diff --git a/setup.py b/setup.py old mode 100644 new mode 100755