ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 5 Pull Requests Actions Releases 14 Activity

Intermediate commit. Ready for some serious testing.

Browse Source
This commit is contained in:
Anne de Jong 2022-06-29 12:25:32 +02:00
parent b561c83448
commit 7095f9d5e7
46 changed files with 1182 additions and 2193 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -32,3 +32,4 @@ lasp/device/lasp_daq.cxx
lasp/c/lasp_config.h lasp/c/lasp_config.h
compile_commands.json compile_commands.json
.cache .cache
lasp_config.h

3
.gitmodules vendored
View File

@ -7,3 +7,6 @@
[submodule "DebugTrace-cpp"] [submodule "DebugTrace-cpp"]
path = DebugTrace-cpp path = DebugTrace-cpp
url = https://github.com/MasatoKokubo/DebugTrace-cpp url = https://github.com/MasatoKokubo/DebugTrace-cpp
[submodule "armadillo-code"]
path = armadillo-code
url = https://gitlab.com/conradsnicta/armadillo-code

59
CMakeLists.txt
View File

@ -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 # To allow linking to static libs from other directories
cmake_policy(SET CMP0079 NEW) cmake_policy(SET CMP0079 NEW)
@ -6,28 +6,25 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/lasp/cma
include("BuildType") include("BuildType")
# This is used for code completion in vim # 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) project(LASP LANGUAGES C CXX)
set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED)
set(CMAKE_C_STANDARD 11) 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_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_RTAUDIO "Compile with RtAudio Daq backend" ON)
option(LASP_HAS_ULDAQ "Compile with UlDaq 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_TRACERNAME "defaulttracer" CACHE STRING "Name of tracer variable containing level")
set(LASP_FFT_BACKEND "FFTW" CACHE STRING "FFT Library backend") set(LASP_FFT_BACKEND "FFTW" CACHE STRING "FFT Library backend")
set_property(CACHE LASP_FFT_BACKEND PROPERTY STRINGS "FFTW" "FFTPack" "None") 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) find_package(pybind11 CONFIG REQUIRED)
# Required for PYBIND11 # Required for PYBIND11
set(POSITION_INDEPENDENT_CODE True) set(POSITION_INDEPENDENT_CODE True)
@ -44,11 +41,12 @@ else()
set(LASP_THREADING_LIBRARIES "") set(LASP_THREADING_LIBRARIES "")
endif() endif()
if(LASP_USE_BLAS)
# link openblas # link openblas
set(BLA_VENDOR OpenBLAS) set(BLA_VENDOR OpenBLAS)
find_package(BLAS REQUIRED) find_package(BLAS REQUIRED)
endif()
# Armadillo
include_directories(SYSTEM armadillo-code/include)
# Reasonable maximum to the nfft size, at 48kHz this is 700s of data... # Reasonable maximum to the nfft size, at 48kHz this is 700s of data...
@ -67,6 +65,7 @@ endif()
# General make flags # General make flags
set(CMAKE_POSITION_INDEPENDENT_CODE ON) set(CMAKE_POSITION_INDEPENDENT_CODE ON)
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wextra -Wno-type-limits \ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wextra -Wno-type-limits \
-Werror=implicit-function-declaration -Wno-unused-parameter \ -Werror=implicit-function-declaration -Wno-unused-parameter \
-Werror=return-type -Wfatal-errors") -Werror=return-type -Wfatal-errors")
@ -98,7 +97,6 @@ else() # Linux compile
include_directories(/usr/local/include/rtaudio) include_directories(/usr/local/include/rtaudio)
include_directories(/usr/include/rtaudio) include_directories(/usr/include/rtaudio)
link_directories(/usr/local/lib) link_directories(/usr/local/lib)
# This should become optional later on, and be added to the windows list as # This should become optional later on, and be added to the windows list as
# well. # well.
@ -108,21 +106,12 @@ endif()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-type-limits") 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 \ set(CMAKE_C_FLAGS_RELEASE "-O2 -mfpmath=sse -march=x86-64 -mtune=native \
-fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions") -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions")
# ############################# End compilation flags # ############################# 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 # Add FFTpack dir if used as FFT backend
if(LASP_FFTPACK_BACKEND) if(LASP_FFTPACK_BACKEND)
add_subdirectory(fftpack) add_subdirectory(fftpack)
@ -131,10 +120,10 @@ if(LASP_FFTPACK_BACKEND)
) )
endif() endif()
include_directories(lasp/c)
include_directories(STL-Threadsafe/include) include_directories(STL-Threadsafe/include)
include_directories(gsl-lite/include) include_directories(gsl-lite/include)
include_directories(DebugTrace-cpp/include) include_directories(DebugTrace-cpp/include)
add_subdirectory(lasp) add_subdirectory(lasp)
add_subdirectory(gsl-lite) add_subdirectory(gsl-lite)
add_subdirectory(test) add_subdirectory(test)
@ -148,21 +137,21 @@ set(DEPS "${CMAKE_CURRENT_SOURCE_DIR}/*.py"
# "lasp_device_lib") # "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}) # # configure_file(${SETUP_PY_IN} ${SETUP_PY})
add_custom_command(OUTPUT ${OUTPUT} # add_custom_command(OUTPUT ${OUTPUT}
COMMAND ${PYTHON_EXECUTABLE} ${SETUP_PY} build # COMMAND ${PYTHON_EXECUTABLE} ${SETUP_PY} build
COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT} # COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
DEPENDS ${DEPS}) # DEPENDS ${DEPS})
add_custom_target(target ALL DEPENDS ${OUTPUT}) # add_custom_target(target ALL DEPENDS ${OUTPUT})
if(DEFINED INSTALL_DEBUG) # if(DEFINED INSTALL_DEBUG)
set(EXTRA_SETUP_ARG --user -e) # set(EXTRA_SETUP_ARG --user -e)
else() # else()
set(EXTRA_SETUP_ARG "") # set(EXTRA_SETUP_ARG "")
endif() # 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} .)")

1
armadillo-code Submodule

@ -0,0 +1 @@
Subproject commit b5727433d342afca53aca0ee16ecf1caaa661821

13
lasp/CMakeLists.txt
View File

@ -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(c)
add_subdirectory(device) add_subdirectory(device)
add_subdirectory(dsp)
# 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)

1
lasp/__init__.py
View File

@ -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_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 .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_atomic import * # Atomic
from .lasp_imptube import * # TwoMicImpedanceTube from .lasp_imptube import * # TwoMicImpedanceTube
from .lasp_measurement import * # Measurement, scaleBlockSens from .lasp_measurement import * # Measurement, scaleBlockSens

45
lasp/c/CMakeLists.txt
View File

@ -1,30 +1,29 @@
configure_file(lasp_config.h.in lasp_config.h)
if(!LASP_DEBUG) if(!LASP_DEBUG)
SET_SOURCE_FILES_PROPERTIES(lasp_sosfilterbank.c PROPERTIES COMPILE_FLAGS -O3) 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_slm.c PROPERTIES COMPILE_FLAGS -O3)
SET_SOURCE_FILES_PROPERTIES(lasp_eq.c PROPERTIES COMPILE_FLAGS -O3) SET_SOURCE_FILES_PROPERTIES(lasp_eq.c PROPERTIES COMPILE_FLAGS -O3)
endif(!LASP_DEBUG) 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})

6
lasp/c/lasp_fft.h
View File

@ -10,6 +10,9 @@
#define LASP_FFT_H #define LASP_FFT_H
#include "lasp_types.h" #include "lasp_types.h"
#include "lasp_mat.h" #include "lasp_mat.h"
#ifdef __cplusplus
extern "C" {
#endif
/** /**
* Load wisdom data from a NULL-terminated character string. Note that at this * 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); void Fft_free(Fft* fft);
#ifdef __cplusplus
}
#endif
#endif // LASP_FFT_H #endif // LASP_FFT_H
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////

41
lasp/c/lasp_math_raw.h
View File

@ -8,6 +8,10 @@
#pragma once #pragma once
#ifndef LASP_MATH_RAW_H #ifndef LASP_MATH_RAW_H
#define LASP_MATH_RAW_H #define LASP_MATH_RAW_H
#ifdef __cplusplus
extern "C" {
#endif
#include "lasp_config.h" #include "lasp_config.h"
#include "lasp_assert.h" #include "lasp_assert.h"
#include "lasp_tracer.h" #include "lasp_tracer.h"
@ -23,40 +27,6 @@
#error "LASP_USE_BLAS should be set to either 0 or 1" #error "LASP_USE_BLAS should be set to either 0 or 1"
#endif #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 /// Positive infinite
#define d_inf (INFINITY) #define d_inf (INFINITY)
@ -502,6 +472,9 @@ static inline void c_conj_inplace(c res[],us size) {
} }
#endif // LASP_USE_BLAS #endif // LASP_USE_BLAS
} }
#ifdef __cplusplus
}
#endif
#endif // LASP_MATH_RAW_H #endif // LASP_MATH_RAW_H
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////

9
lasp/c/lasp_nprocs.h
View File

@ -3,16 +3,23 @@
// Author: J.A. de Jong - ASCEE // Author: J.A. de Jong - ASCEE
// //
// Description: Implemententation of a function to determine the number // Description: Implemententation of a function to determine the number
// of processors. // of processors. Useful for spreading out computational load over multiple
// cores.
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
#pragma once #pragma once
#ifndef LASP_NPROCS_H #ifndef LASP_NPROCS_H
#define LASP_NPROCS_H #define LASP_NPROCS_H
#include "lasp_types.h" #include "lasp_types.h"
#ifdef __cplusplus
extern C {
#endif
/** /**
* @return The number of SMP processors * @return The number of SMP processors
*/ */
us getNumberOfProcs(); us getNumberOfProcs();
#ifdef __cplusplus
}
#endif
#endif // LASP_NPROCS_H #endif // LASP_NPROCS_H

472
lasp/c/lasp_siggen.c
View File

@ -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<N; n++) {
setvecval(data, n, d_sin(phase));
d fn = fl + ((d) n)/N*(fu + eps -fl);
phase += 2*number_pi*Dt*fn;
}
}
else {
/* Continous sweep */
TRACE(15, "continuous sweep");
iVARTRACE(17, N);
dVARTRACE(17, fl);
dVARTRACE(17, fu);
const us Nf = N/2;
const us Nb = N-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);
for(us n = 0; n<=N; n++) {
/* iVARTRACE(17, n); */
if(n<N) {
setvecval(data, 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;
/* 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<N; n++) {
setvecval(data, n, d_sin(phase));
d fn = fl*d_pow(k,((d) n)/N);
phase += 2*number_pi*Dt*fn;
}
} else {
TRACE(15, "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;
dVARTRACE(15, 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;
dVARTRACE(15, k);
}
iVARTRACE(15, K);
dVARTRACE(15, k1);
dVARTRACE(15, k);
dVARTRACE(15, E);
for(us n = 0; n<=N; n++) {
/* iVARTRACE(17, n); */
if(n<N) {
setvecval(data, 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); */
/* 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; i<samples->n_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);
}
//////////////////////////////////////////////////////////////////////

100
lasp/c/lasp_siggen.h
View File

@ -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
//////////////////////////////////////////////////////////////////////

78
lasp/c/lasp_tracer.c
View File

@ -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 <stdio.h>
#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<ASCEE_FN_LEVEL;i++) {
printf("--");
}
printf("* ");
}
void trace_impl(const char* file,int pos, const char * string){
indent_trace();
printf(annestr(TRACERNAME) ":%s:%i: %s\n",file,pos,string);
}
void fstrace_impl(const char* file,int pos,const char* fn){
ASCEE_FN_LEVEL++;
indent_trace();
printf(annestr(TRACERNAME) ":%s:%i: start function: %s()\n",file,pos,fn);
}
void fetrace_impl(const char* file,int pos,const char* fn){
indent_trace();
printf(annestr(TRACERNAME) ":%s:%i: end function: %s()\n",file,pos,fn);
ASCEE_FN_LEVEL--;
}
void ivartrace_impl(const char* pos,int line,const char* varname, int var){
indent_trace();
printf(annestr(TRACERNAME) ":%s:%i: %s = %i\n",pos,line,varname,var);
}
void uvartrace_impl(const char* pos,int line,const char* varname,size_t var){
indent_trace();
printf(annestr(TRACERNAME) ":%s:%i: %s = %zu\n",pos,line,varname,var);
}
void dvartrace_impl(const char* pos,int line,const char* varname, d var){
indent_trace();
printf(annestr(TRACERNAME) ":%s:%i: %s = %0.5e\n",pos,line,varname,var);
}
void cvartrace_impl(const char* pos,int line,const char* varname, c var){
indent_trace();
printf(annestr(TRACERNAME) ":%s:%i: %s = %0.5e+%0.5ei\n",pos,line,varname,creal(var),cimag(var));
}
#endif
//////////////////////////////////////////////////////////////////////

238
lasp/c/lasp_tracer.h
View File

@ -1,238 +0,0 @@
// ascee_tracer.h
//
// Author: J.A. de Jong - ASCEE
//
// Description:
// Basic tracing code for debugging.
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_TRACER_H
#define LASP_TRACER_H
#include "lasp_config.h"
#include "lasp_types.h"
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#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
//////////////////////////////////////////////////////////////////////

111
lasp/c/lasp_window.c
View File

@ -1,111 +0,0 @@
// window.c
//
// Author: J.A. de Jong -ASCEE
//
// Description:
//
//////////////////////////////////////////////////////////////////////
#include "lasp_window.h"
#include "lasp_signals.h"
#include <stdlib.h>
/**
* Compute the Hann window
*
* @param i index
* @param N Number of indices
*/
static d hann(us i,us N) {
dbgassert(i<N,"Invalid index for window function Hann");
return d_pow(d_sin(number_pi*i/(N-1)),2);
}
/**
* Compute the Hamming window
*
* @param i index
* @param N Number of indices
*/
static d hamming(us i,us N) {
dbgassert(i<N,"Invalid index for window function Hamming");
d alpha = 25.0/46.0;
return alpha-(1-alpha)*d_cos(2*number_pi*i/(N-1));
}
/**
* Compute the Blackman window
*
* @param i index
* @param N Number of indices
*/
static d blackman(us i,us N) {
dbgassert(i<N,"Invalid index for window function Blackman");
d a0 = 7938./18608.;
d a1 = 9240./18608.;
d a2 = 1430./18608.;
return a0-a1*d_cos(2*number_pi*i/(N-1))+a2*d_cos(4*number_pi*i/(N-1));
}
/**
* Compute the Rectangular window
*
* @param i index
* @param N Number of indices
*/
static d rectangle(us i,us N) {
dbgassert(i<N,"Invalid index for window function Hann");
return 1.0;
}
static d bartlett(us n,us N) {
dbgassert(n<N,"Invalid index for window function Bartlett");
return 1 - d_abs(2*(n - (N-1)/2.)/(N-1));
}
int window_create(const WindowType wintype,vd* result,d* win_pow) {
fsTRACE(15);
dbgassert(result && win_pow,NULLPTRDEREF);
assert_vx(result);
us nfft = result->n_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<nfft;index++) {
/* Compute the window function value */
d val = win_fun(index,nfft);
/* Set the value in the vector */
setvecval(result,index,val);
}
/* Store window power in result */
*win_pow = signal_power(result);
feTRACE(15);
return LASP_SUCCESS;
}
//////////////////////////////////////////////////////////////////////

37
lasp/c/lasp_window.h
View File

@ -1,37 +0,0 @@
// window.h
//
// Author: J.A. de Jong - ASCEE
//
// Description:
//
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_WINDOW_H
#define LASP_WINDOW_H
#include "lasp_mat.h"
typedef enum {
Hann = 0,
Hamming = 1,
Rectangular = 2,
Bartlett = 3,
Blackman = 4,
} WindowType;
/**
* Create a Window function, store it in the result
*
* @param[in] wintype Enumerated type, corresponding to the window
* function.
* @param[out] result Vector where the window values are stored
* @param[out] win_power Here, the overall power of the window will be
* returned.
*
* @return status code, SUCCESS on success.
*/
int window_create(const WindowType wintype,vd* result,d* win_power);
#endif // LASP_WINDOW_H
//////////////////////////////////////////////////////////////////////

10
lasp/config.pxi.in
View File

@ -1,10 +0,0 @@
import numpy as np
cimport numpy as cnp
# Do this, our segfaults will be your destination
cnp.import_array()
DEF LASP_DOUBLE_PRECISION = "@LASP_DOUBLE_PRECISION@"
DEF LASP_DEBUG_CYTHON = "@LASP_DEBUG@"
from libcpp cimport bool

30
lasp/device/CMakeLists.txt
View File

@ -1,13 +1,23 @@
file(GLOB device_files *.cpp) add_library(lasp_device_lib OBJECT
# set(cpp_daq_linklibs ${LASP_THREADING_LIBRARIES}) lasp_daq.cpp
lasp_daqconfig.cpp
lasp_daqdata.cpp
lasp_deviceinfo.cpp
lasp_rtaudiodaq.cpp
lasp_streammgr.cpp
lasp_uldaq.cpp
)
# if(win32) target_include_directories(lasp_device_lib PUBLIC ../dsp)
# list(APPEND lasp_device_linklibs python${python_version_windll}) target_include_directories(lasp_device_lib PUBLIC ../c)
# endif(win32) if(LASP_HAS_ULDAQ)
target_link_libraries(lasp_device_lib uldaq)
endif()
if(LASP_HAS_RTAUDIO)
target_link_libraries(lasp_device_lib rtaudio)
endif()
target_link_libraries(lasp_device_lib lasp_dsp_lib)
add_library(lasp_device_lib OBJECT ${device_files}) pybind11_add_module(lasp_device lasp_devicepybind.cpp)
target_link_libraries(lasp_device PRIVATE lasp_device_lib)
# pybind11_add_module(lasp_device lasp_pydevicepybid.cpp)
# target_link_libraries(lasp_device PRIVATE lasp_device_lib)

5
lasp/device/__init__.py
View File

@ -1,5 +1,2 @@
from .lasp_device_common import * from .lasp_device import *
from .lasp_deviceinfo import *
from .lasp_daqconfig import *
from .lasp_daq import *

3
lasp/device/lasp_daq.cpp
View File

@ -70,6 +70,9 @@ double Daq::samplerate() const {
DataTypeDescriptor::DataType Daq::dataType() const { DataTypeDescriptor::DataType Daq::dataType() const {
return availableDataTypes.at(dataTypeIndex); return availableDataTypes.at(dataTypeIndex);
} }
DataTypeDescriptor Daq::dtypeDescr() const {
return dtype_map.at(dataType());
}
double Daq::inputRangeForChannel(us ch) const { double Daq::inputRangeForChannel(us ch) const {
if (!(ch < ninchannels)) { if (!(ch < ninchannels)) {

32
lasp/device/lasp_daq.h
View File

@ -1,11 +1,29 @@
#pragma once #pragma once
#include <memory>
#include "lasp_config.h" #include "lasp_config.h"
#include "lasp_daqdata.h" #include "lasp_daqdata.h"
#include "lasp_deviceinfo.h" #include "lasp_deviceinfo.h"
#include "lasp_types.h" #include "lasp_types.h"
#include <functional> #include <functional>
#include <gsl/gsl-lite.hpp> #include <gsl/gsl-lite.hpp>
#include <memory>
/**
* @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 * @brief Callback of DAQ for input data. Callback should return
@ -48,7 +66,7 @@ public:
* data is presented, the function is called with a nullptr as argument. * data is presented, the function is called with a nullptr as argument.
*/ */
virtual void start(std::optional<InDaqCallback> inCallback, virtual void start(std::optional<InDaqCallback> inCallback,
std::optional<OutDaqCallback> outCallback); std::optional<OutDaqCallback> outCallback) = 0;
/** /**
* @brief Stop the Daq device. Throws an exception if the device is not * @brief Stop the Daq device. Throws an exception if the device is not
@ -63,7 +81,7 @@ public:
*/ */
virtual bool isRunning() const = 0; virtual bool isRunning() const = 0;
virtual ~Daq(){}; virtual ~Daq() = default;
/** /**
* @brief Returns the number of enabled input channels * @brief Returns the number of enabled input channels
* *
@ -118,6 +136,14 @@ public:
us framesPerBlock() const { us framesPerBlock() const {
return availableFramesPerBlock.at(framesPerBlockIndex); 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 * @brief Whether the device runs in duplex mode (both input and output), or
* false if only input / only output. * false if only input / only output.

6
lasp/device/lasp_daqconfig.h
View File

@ -1,8 +1,8 @@
#pragma once #pragma once
#include "lasp_config.h"
#include "lasp_types.h"
#include <map> #include <map>
#include <vector> #include <vector>
#include "lasp_config.h"
#include "lasp_types.h"
using std::string; using std::string;
using std::to_string; using std::to_string;
@ -200,7 +200,7 @@ public:
* @param deviceinfo DeviceInfo structure * @param deviceinfo DeviceInfo structure
*/ */
DaqConfiguration(const DeviceInfo &DeviceInfo); DaqConfiguration(const DeviceInfo &DeviceInfo);
DaqConfiguration() = delete; DaqConfiguration();
/** /**
* @brief Check to see whether the DAQ configuration matches with the device. * @brief Check to see whether the DAQ configuration matches with the device.

4
lasp/device/lasp_daqdata.h
View File

@ -1,10 +1,10 @@
#pragma once #pragma once
#include "lasp_daqconfig.h"
#include "lasp_types.h"
#include <functional> #include <functional>
#include <gsl/gsl-lite.hpp> #include <gsl/gsl-lite.hpp>
#include <memory> #include <memory>
#include <optional> #include <optional>
#include "lasp_daqconfig.h"
#include "lasp_types.h"
/** /**
* @brief Data coming from / going to DAQ. **Non-interleaved format**, which * @brief Data coming from / going to DAQ. **Non-interleaved format**, which

161
lasp/device/lasp_devicepybind.cpp
View File

@ -1,31 +1,33 @@
#if 0
#include "lasp_devicepybind.h"
#include <algorithm> #include <algorithm>
#include <assert.h> #include <assert.h>
#include <numpy/arrayobject.h>
#include <pybind11/numpy.h> #include <pybind11/numpy.h>
#include <stdint.h> #include <stdint.h>
#include <thread> #include <iostream>
#include "lasp_daqconfig.h"
#include "lasp_deviceinfo.h"
using std::cerr; using std::cerr;
namespace py = pybind11; namespace py = pybind11;
PYBIND11_MODULE(lasp_daq, m) { PYBIND11_MODULE(lasp_device, m) {
m.doc() = "Lasp DAQ interface"; m.doc() = "Lasp device interface";
/// DataType /// DataType
py::class_<DataType> datatype(m, "DataType"); py::class_<DataTypeDescriptor> dtype_desc(m, "DataTypeDescriptor");
datatype.def_readonly("name", &DataType::name); dtype_desc.def_readonly("name", &DataTypeDescriptor::name);
datatype.def_readonly("sw", &DataType::sw); dtype_desc.def_readonly("sw", &DataTypeDescriptor::sw);
datatype.def_readonly("is_floating", &DataType::is_floating); dtype_desc.def_readonly("is_floating", &DataTypeDescriptor::is_floating);
py::enum_<DataTypeDescriptor::DataType>(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 /// DaqApi
py::class_<DaqApi> daqapi(m, "DaqApi"); py::class_<DaqApi> daqapi(m, "DaqApi");
@ -105,132 +107,5 @@ PYBIND11_MODULE(lasp_daq, m) {
daqconfig.def("match", &DaqConfiguration::match); daqconfig.def("match", &DaqConfiguration::match);
/// Daq
py::class_<PyDaq> daq(m, "Daq");
daq.def(py::init<DeviceInfo,DaqConfiguration>());
/* 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<float>::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<void *> *inQueue = input ? &_inQueue : nullptr;
SafeQueue<void *> *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<d>(
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

44
lasp/device/lasp_devicepybind.h
View File

@ -1,44 +0,0 @@
#pragma once
#include "lasp_daq.h"
#include <atomic>
#include <vector>
#include <pybind11/pybind11.h>
namespace py = pybind11;
namespace std {
class thread;
};
class PyDaq {
Daq *_daq = nullptr;
std::thread *_daqThread = nullptr;
SafeQueue<void *> _inQueue;
SafeQueue<void *> _outQueue;
py::function _audioCallback;
std::atomic<bool> _stopThread;
std::atomic<bool> _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<DeviceInfo> getDeviceInfo();
const Daq &daq() const { return *_daq; }
private:
void threadFcn();
};

12
lasp/device/lasp_rtaudiodaq.cpp
View File

@ -300,14 +300,10 @@ public:
return 0; return 0;
} }
~RtAudioDaq() { // RtAudio documentation says: if a stream is open, it will be stopped and
if (isRunning()) { // closed automatically on deletion. Therefore the destructor here is a
stop(); // default one.
} ~RtAudioDaq() = default;
if (rtaudio.isStreamOpen()) {
rtaudio.closeStream();
}
}
}; };
std::unique_ptr<Daq> createRtAudioDevice(const DeviceInfo &devinfo, std::unique_ptr<Daq> createRtAudioDevice(const DeviceInfo &devinfo,

151
lasp/device/lasp_streammgr.cpp Normal file
View File

@ -0,0 +1,151 @@
#include <assert.h>
#include <algorithm>
#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<std::mutex> 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 <typename T> bool fillData(DaqData &data, const vd &signal) {
assert(data.nframes == signal.size());
T* res = reinterpret_cast<T*>(data.raw_ptr());
if (std::is_floating_point<T>()) {
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<T>::max());
}
}
}
return true;
}
void StreamMgr::setSiggen(std::shared_ptr<Siggen> siggen) {
std::scoped_lock<std::mutex> 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<std::mutex> lck(_siggen_mtx);
if(_siggen) {
vd signal = _siggen->genSignal(data.nframes);
switch (data.dtype) {
case (DataTypeDescriptor::DataType::dtype_fl32):
fillData<float>(data, signal);
break;
case (DataTypeDescriptor::DataType::dtype_fl64):
fillData<double>(data, signal);
break;
case (DataTypeDescriptor::DataType::dtype_int8):
fillData<int8_t>(data, signal);
break;
case (DataTypeDescriptor::DataType::dtype_int16):
fillData<int16_t>(data, signal);
break;
case (DataTypeDescriptor::DataType::dtype_int32):
fillData<int32_t>(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<std::mutex> lck(_inDataHandler_mtx);
_inDataHandlers.push_back(&handler);
}
void StreamMgr::removeInDataHandler(InDataHandler &handler) {
std::scoped_lock<std::mutex> lck(_inDataHandler_mtx);
_inDataHandlers.remove(&handler);
}

187
lasp/device/lasp_streammgr.h Normal file
View File

@ -0,0 +1,187 @@
#pragma once
#include "lasp_daq.h"
#include "lasp_siggen.h"
#include <list>
#include <memory>
#include <mutex>
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<StreamType, std::unique_ptr<Daq>> _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<InDataHandler *> _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> _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<Siggen> s);
};

1
lasp/device/lasp_uldaq.cpp
View File

@ -8,7 +8,6 @@
#include <iostream> #include <iostream>
#include <memory> #include <memory>
#include <mutex> #include <mutex>
#include <span>
#include <stdexcept> #include <stdexcept>
#include <thread> #include <thread>
#include <uldaq.h> #include <uldaq.h>

21
lasp/dsp/CMakeLists.txt Normal file
View File

@ -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})

1
lasp/dsp/__init__.py Normal file
View File

@ -0,0 +1 @@
from .lasp_dsp import *

0
lasp/c/lasp_config.h.in → lasp/dsp/lasp_config.h.in
View File

20
lasp/dsp/lasp_dsp_pybind.cpp Normal file
View File

@ -0,0 +1,20 @@
#include "lasp_window.h"
#include <iostream>
#include <pybind11/pybind11.h>
using std::cerr;
namespace py = pybind11;
PYBIND11_MODULE(lasp_dsp, m) {
py::class_<Window> w(m, "Window");
py::enum_<Window::WindowType>(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();
}

0
lasp/dsp/lasp_filter.cpp Normal file
View File

14
lasp/dsp/lasp_filter.h Normal file
View File

@ -0,0 +1,14 @@
#pragma once
#include <vector>
#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;
};

45
lasp/dsp/lasp_mathtypes.h Normal file
View File

@ -0,0 +1,45 @@
#pragma once
#include <armadillo>
#include "lasp_types.h"
#include <cmath>
#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<d>;
using dmat = arma::Mat<d>;
const d number_pi = arma::datum::pi;

42
lasp/dsp/lasp_siggen.cpp Normal file
View File

@ -0,0 +1,42 @@
#include "lasp_siggen.h"
#include "lasp_mathtypes.h"
#include <cassert>
#include <type_traits>
inline d level_amp(d level_dB){
return pow(10, level_dB/20);
}
using mutexlock = std::scoped_lock<std::mutex>;
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>& 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();
}

77
lasp/dsp/lasp_siggen.h Normal file
View File

@ -0,0 +1,77 @@
#pragma once
#include <memory>
#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> _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<Filter>& 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);
};

245
lasp/dsp/lasp_siggen_impl.cpp Normal file
View File

@ -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<vd>(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);
}
}
}

79
lasp/dsp/lasp_siggen_impl.h Normal file
View File

@ -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;
};

8
lasp/c/lasp_types.h → lasp/dsp/lasp_types.h
View File

@ -51,14 +51,6 @@ typedef double complex c;
#endif #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 #endif // LASP_TYPES_H
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////

57
lasp/dsp/lasp_window.cpp Normal file
View File

@ -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;
}
}

73
lasp/dsp/lasp_window.h Normal file
View File

@ -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);
};

12
lasp/lasp_common.py
View File

@ -2,11 +2,11 @@
import shelve, logging, sys, appdirs, os, platform import shelve, logging, sys, appdirs, os, platform
import numpy as np import numpy as np
from .wrappers import Window as wWindow
from collections import namedtuple from collections import namedtuple
from dataclasses import dataclass from dataclasses import dataclass
from dataclasses_json import dataclass_json from dataclasses_json import dataclass_json
from enum import Enum, unique, auto from enum import Enum, unique, auto
from .dsp import Window as wWindow
""" """
Common definitions used throughout the code. Common definitions used throughout the code.
@ -271,11 +271,11 @@ class this_lasp_shelve(Shelve):
@unique @unique
class Window(Enum): class Window(Enum):
hann = (wWindow.hann, 'Hann') hann = (wWindow.Hann, 'Hann')
hamming = (wWindow.hamming, 'Hamming') hamming = (wWindow.Hamming, 'Hamming')
rectangular = (wWindow.rectangular, 'Rectangular') rectangular = (wWindow.Rectangular, 'Rectangular')
bartlett = (wWindow.bartlett, 'Bartlett') bartlett = (wWindow.Bartlett, 'Bartlett')
blackman = (wWindow.blackman, 'Blackman') blackman = (wWindow.Blackman, 'Blackman')
@staticmethod @staticmethod
def fillComboBox(cb): def fillComboBox(cb):

798
lasp/wrappers.pyx
View File

@ -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(<const char*> &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,<WindowType> 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,
<WindowType> 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)

0
setup.py Normal file → Executable file
View File