Major cleanup and improvements of filter design code.

.gitignore

@@ -1,5 +1,6 @@
 *.a
 *.pyc
+lasp_rtaudio.cxx
 dist
 CMakeFiles
 CMakeCache.txt
@@ -23,6 +24,5 @@ LASP.egg-info
 lasp_octave_fir.*
 lasp/ui_*
 resources_rc.py
-lasp/device/lasp_daqdevice.c
 .ropeproject
 .spyproject

CMakeLists.txt

@@ -1,5 +1,8 @@
-cmake_minimum_required (VERSION 3.0)
-project(beamforming)
+cmake_minimum_required (VERSION 3.0) 
+
+# This is used for code completion in vim
+set(CMAKE_EXPORT_COMPILE_COMMANDS=ON)
+project(LASP)
 
 # Whether we want to use blas yes or no
 set(LASP_USE_BLAS TRUE)
@@ -20,17 +23,32 @@ add_definitions(-DLASP_MAX_NFFT=33554432) # 2**25
 
 add_definitions(-D_REENTRANT)
 
+# ############### Choose an fft backend here
+#set(LASP_FFT_BACKEND fftpack)
+set(LASP_FFT_BACKEND "fftw")
+
+if(LASP_FFT_BACKEND STREQUAL "fftw")
+    find_library(FFTW_LIBRARY NAMES fftw3 fftw)
+    set(FFTW_LIBRARIES "${FFTW_LIBRARY}")
+    set(LASP_FFT_LIBRARY "${FFTW_LIBRARIES}")
+    add_definitions(-DLASP_FFT_BACKEND_FFTW)
+elseif(LASP_FFT_BACKEND STREQUAL "fftpack")
+    add_definitions(-DLASP_FFT_BACKEND_FFTPACK)
+    set(LASP_FFT_LIBRARY "fftpack")
+endif()
+
 
 if(LASP_FLOAT STREQUAL "double")
-  add_definitions(-DLASP_FLOAT=64)
-  add_definitions(-DLASP_DOUBLE_PRECISION)
+    add_definitions(-DLASP_FLOAT=64)
+    add_definitions(-DLASP_DOUBLE_PRECISION)
 else()
-  add_definitions(-DLASP_FLOAT=32)
+    add_definitions(-DLASP_FLOAT=32)
+    add_definitions(-DLASP_SINGLE_PRECISION)
 endif(LASP_FLOAT STREQUAL "double")
 
 if(NOT DEFINED LASP_DEBUG)
-  message(FATAL_ERROR "LASP_DEBUG flag not defined. Please set -DLASP_DEBUG=TRUE
-   or -DLASP_DEBUG=FALSE")
+    message(FATAL_ERROR "LASP_DEBUG flag not defined. Please set -DLASP_DEBUG=TRUE
+    or -DLASP_DEBUG=FALSE")
 endif(NOT DEFINED LASP_DEBUG)
 
 # ##################### END Cmake variables converted to a macro
@@ -38,64 +56,62 @@ set(Python_ADDITIONAL_VERSIONS "3")
 # #################### Setting definitions and debug-specific compilation flags
 
 if(LASP_DEBUG)
-  set(TRACERNAME LASPTracer)
-  set(CMAKE_BUILD_TYPE Debug)
-  message("Building debug code")
-  set(CMAKE_BUILD_TYPE Debug)
-  add_definitions(-DLASP_DEBUG=1)
-  add_definitions(-DTRACERNAME=${TRACERNAME})
-  add_definitions(-DDEBUG)
-  add_definitions(-DTRACER=1)
-
-  set(CYTHON_VARIABLES "#cython: boundscheck=True")
-  # This will produce html files
-  set(CYTHON_ANNOTATE ON)
-  # Add the __FILENAME__ macro
-  # set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -D__FILENAME__='\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"'")
+    set(TRACERNAME LASPTracer)
+    set(LASP_DEBUG_CYTHON=True)
+    set(CMAKE_BUILD_TYPE Debug)
+    message("Building debug code")
+    set(CMAKE_BUILD_TYPE Debug)
+    add_definitions(-DLASP_DEBUG=1)
+    add_definitions(-DTRACERNAME=${TRACERNAME})
+    add_definitions(-DDEBUG)
+    add_definitions(-DTRACER=1)
+    # This will produce html files
+    set(CYTHON_ANNOTATE ON)
+    # Add the __FILENAME__ macro
+    # set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -D__FILENAME__='\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"'")
 else()
-  message("Building LASP for release")
-  set(CMAKE_BUILD_TYPE Release)
-  set(CYTHON_ANNOTATE OFF)
-  set(CYTHON_NO_DOCSTRINGS ON)
-  set(CYTHON_VARIABLES "#cython: boundscheck=False,wraparound=False,embedsignature=False,emit_code_comments=False")
-  # Strip unnecessary symbols
-  # set(CMAKE_SHARED_LINKER_FLAGS "-Wl,--gc-sections")
-  # set(CMAKE_MODULE_LINKER_FLAGS "-Wl,--gc-sections")
+    message("Building LASP for release")
+    set(CMAKE_BUILD_TYPE Release)
+    set(LASP_DEBUG_CYTHON=False)
+    set(CYTHON_ANNOTATE OFF)
+    set(CYTHON_NO_DOCSTRINGS ON)
+    # Strip unnecessary symbols
+    # set(CMAKE_SHARED_LINKER_FLAGS "-Wl,--gc-sections")
+    # set(CMAKE_MODULE_LINKER_FLAGS "-Wl,--gc-sections")
 
-  add_definitions(-DTRACER=0 -DNDEBUG)
+    add_definitions(-DTRACER=0 -DNDEBUG)
 endif(LASP_DEBUG)
 
 # The last argument here takes care of calling SIGABRT when an integer overflow
 # occures.
 ############################## General compilation flags (independent of debug mode, windows or linux)
-set(CYTHON_EXTRA_C_FLAGS "-Wno-sign-compare -Wno-cpp -Wno-implicit-fallthrough\
- -Wno-incompatible-pointer-types -Wno-strict-aliasing")
+set(CYTHON_EXTRA_C_FLAGS "-Wno-sign-compare -Wno-cpp -Wno-implicit-fallthrough -Wno-incompatible-pointer-types -Wno-strict-aliasing")
 
 # General make flags
 set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fPIC -std=c11  -Wall -Wextra -Wno-type-limits \
- -Werror=implicit-function-declaration -Werror=incompatible-pointer-types \
+-Werror=implicit-function-declaration -Werror=incompatible-pointer-types \
 -Werror=return-type")
 
 # Debug make flags
 set(CMAKE_C_FLAGS_DEBUG "-g" )
 
 set(CMAKE_C_FLAGS_RELEASE  "-O2 -mfpmath=sse -march=x86-64 -mtune=native \
-  -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions")
+-fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions")
 
 # set(CMAKE_C_FLAGS_RELEASE  "-O2 -march=native -mtune=native -fomit-frame-pointer")
 
 if(LASP_USE_BLAS)
-  add_definitions(-DLASP_USE_BLAS=1)
+    add_definitions(-DLASP_USE_BLAS=1)
 else()
-  add_definitions(-DLASP_USE_BLAS=0)
+    add_definitions(-DLASP_USE_BLAS=0)
 endif(LASP_USE_BLAS)
 
 # ############################# End compilation flags
 
 if(CMAKE_SYSTEM_NAME STREQUAL "Windows")
-  set(win32 true)
+    set(win32 true)
 else()
-  set(win32 false)
+    set(win32 false)
 endif(CMAKE_SYSTEM_NAME STREQUAL "Windows")
 
 find_package(PythonLibs REQUIRED )
@@ -103,9 +119,9 @@ find_package(PythonInterp REQUIRED)
 
 add_subdirectory(fftpack)
 include_directories(
-  fftpack
-  lasp/c
-  )
+    fftpack
+    lasp/c
+    )
 
 add_subdirectory(lasp)
 add_subdirectory(test)
@@ -114,18 +130,18 @@ add_subdirectory(test)
 set(SETUP_PY    "${CMAKE_CURRENT_BINARY_DIR}/setup.py")
 
 set(DEPS  "${CMAKE_CURRENT_SOURCE_DIR}/*.py"
-          "${CMAKE_CURRENT_SOURCE_DIR}/lasp/*.py"
-          "wrappers"
-          "lasp_daqdevice")
+    "${CMAKE_CURRENT_SOURCE_DIR}/lasp/*.py"
+    "wrappers"
+    "lasp_rtaudio")
 
-             # )
+# )
 set(OUTPUT      "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp")
 
 # configure_file(${SETUP_PY_IN} ${SETUP_PY})
 add_custom_command(OUTPUT ${OUTPUT}
-                   COMMAND ${PYTHON_EXECUTABLE} ${SETUP_PY} build
-                   COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
-                   DEPENDS ${DEPS})
+    COMMAND ${PYTHON_EXECUTABLE} ${SETUP_PY} build
+    COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+    DEPENDS ${DEPS})
 
 add_custom_target(target ALL DEPENDS ${OUTPUT})
 
@@ -133,7 +149,7 @@ if(DEFINED INSTALL_DEBUG)
     set(EXTRA_SETUP_ARG --user -e)
 else()
     set(EXTRA_SETUP_ARG "")
- endif()
+endif()
 
 
 install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pip install ${EXTRA_SETUP_ARG} .)")

lasp/CMakeLists.txt

@@ -1,6 +1,7 @@
 
 configure_file(config.pxi.in config.pxi)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
+set(CYTHON_EXECUTABLE "cython3")
 include(UseCython)
 find_package(Numpy REQUIRED )
 
@@ -9,8 +10,8 @@ include_directories(
    .
    c
  )
- add_subdirectory(c)
- add_subdirectory(device)
+add_subdirectory(c)
+add_subdirectory(device)
 
 set_source_files_properties(wrappers.c PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} ${CYTHON_EXTRA_C_FLAGS}")
 cython_add_module(wrappers wrappers.pyx)

lasp/c/CMakeLists.txt

@@ -1,6 +1,5 @@
 if(!LASP_DEBUG)
-# SET_SOURCE_FILES_PROPERTIES(y.c PROPERTIES COMPILE_FLAGS -O3)
-# SET_SOURCE_FILES_PROPERTIES(x.c PROPERTIES COMPILE_FLAGS -O3)
+SET_SOURCE_FILES_PROPERTIES(lasp_sosfilterbank.c PROPERTIES COMPILE_FLAGS -O3)
 endif(!LASP_DEBUG)
 
 add_library(lasp_lib
@@ -14,14 +13,15 @@ add_library(lasp_lib
   lasp_aps.c
   lasp_ps.c
   lasp_mq.c
+  lasp_siggen.c
   lasp_worker.c
   lasp_dfifo.c
-  lasp_filterbank.c
-  # lasp_octave_fir.c
+  lasp_firfilterbank.c
+  lasp_sosfilterbank.c
   lasp_decimation.c
   lasp_sp_lowpass.c
   )
 
 
 
-target_link_libraries(lasp_lib fftpack openblas)
+target_link_libraries(lasp_lib ${LASP_FFT_LIBRARY} openblas)

lasp/c/lasp_assert.h

@@ -30,6 +30,8 @@ void DBG_AssertFailedExtImplementation(const char* file,
     	DBG_AssertFailedExtImplementation(__FILE__, __LINE__, assert_string );  \
     }
 
+#define assertvalidptr(ptr) dbgassert(ptr,NULLPTRDEREF)
+
 #else // LASP_DEBUG not defined
 
 #define dbgassert(assertion, assert_string)

lasp/c/lasp_decimation.c

@@ -6,7 +6,7 @@
 // Implementation of the decimator
 //////////////////////////////////////////////////////////////////////
 #include "lasp_decimation.h"
-#include "lasp_filterbank.h"
+#include "lasp_firfilterbank.h"
 #include "lasp_tracer.h"
 #include "lasp_alloc.h"
 #include "lasp_dfifo.h"
@@ -34,7 +34,7 @@ static __thread DecFilter DecFilters[] = {
 typedef struct Decimator_s {
     us nchannels;
     us dec_fac;
-    FilterBank** fbs;
+    Firfilterbank** fbs;
     dFifo* output_fifo;
 } Decimator;
 
@@ -61,14 +61,14 @@ Decimator* Decimator_create(us nchannels,DEC_FAC df) {
     
     /* Create the filterbanks */
     Decimator* dec = a_malloc(sizeof(Decimator));
-    dec->fbs = a_malloc(sizeof(FilterBank*)*nchannels);
+    dec->fbs = a_malloc(sizeof(Firfilterbank*)*nchannels);
     dec->nchannels = nchannels;
     dec->dec_fac = filter->dec_fac;
 
     dmat h = dmat_foreign_data(filter->ntaps,1,filter->h,false);
 
     for(us channelno=0;channelno<nchannels;channelno++) {
-        dec->fbs[channelno] = FilterBank_create(&h,DEC_FFT_LEN);
+        dec->fbs[channelno] = Firfilterbank_create(&h,DEC_FFT_LEN);
     }
 
     dmat_free(&h);
@@ -128,10 +128,10 @@ dmat Decimator_decimate(Decimator* dec,const dmat* samples) {
                                          chan);
 
         /* Low-pass filter stuff */
-        dmat filtered_res = FilterBank_filter(dec->fbs[chan],
+        dmat filtered_res = Firfilterbank_filter(dec->fbs[chan],
                                                  &samples_channel);
 
-        dbgassert(filtered_res.n_cols == 1,"Bug in FilterBank");
+        dbgassert(filtered_res.n_cols == 1,"Bug in Firfilterbank");
 
         vd_free(&samples_channel);
 
@@ -148,7 +148,7 @@ dmat Decimator_decimate(Decimator* dec,const dmat* samples) {
                                             1);
 
             dbgassert(filtered_res.n_rows == filtered_col.n_rows,
-                      "Not all FilterBank's have same output number"
+                      "Not all Firfilterbank's have same output number"
                       " of rows!");
 
             dmat_copy_rows(&filtered_col,
@@ -209,7 +209,7 @@ void Decimator_free(Decimator* dec) {
     dFifo_free(dec->output_fifo);
 
     for(us chan=0;chan<dec->nchannels;chan++) {
-        FilterBank_free(dec->fbs[chan]);
+        Firfilterbank_free(dec->fbs[chan]);
     }
 
     a_free(dec->fbs);

lasp/c/lasp_fft.c

@@ -9,12 +9,24 @@
 #include "lasp_tracer.h"
 #include "lasp_fft.h"
 #include "lasp_types.h"
-#include "fftpack.h"
 
+#ifdef LASP_FFT_BACKEND_FFTPACK
+#include "fftpack.h"
 typedef struct Fft_s {
     us nfft;
-    vd fft_work;
-} Fft;
+    vd fft_work; // Storage memory for fftpack
+};
+#elif defined LASP_FFT_BACKEND_FFTW
+#include <fftw3.h>
+typedef struct Fft_s {
+    us nfft;
+    fftw_plan forward_plan;
+    fftw_plan reverse_plan;
+    c* complex_storage;
+    d* real_storage;
+};
+
+#endif
 
 Fft* Fft_create(const us nfft) {
     fsTRACE(15);
@@ -24,72 +36,103 @@ Fft* Fft_create(const us nfft) {
     }
 
     Fft* fft = a_malloc(sizeof(Fft));
-    if(fft==NULL) {
-        WARN("Fft allocation failed");
-        return NULL;
-    }
 
     fft->nfft = nfft;
 
+#ifdef LASP_FFT_BACKEND_FFTPACK
     /* Initialize foreign fft lib */
     fft->fft_work = vd_alloc(2*nfft+15);
-
     npy_rffti(nfft,getvdval(&fft->fft_work,0));
     check_overflow_vx(fft->fft_work);
+#elif defined LASP_FFT_BACKEND_FFTW
+    fft->complex_storage = fftw_malloc(sizeof(c) * (nfft/2 + 1));
+    fft->real_storage = fftw_malloc(sizeof(d) * nfft);
+
+    fft->forward_plan = fftw_plan_dft_r2c_1d(nfft,
+            fft->real_storage,
+            fft->complex_storage,
+            FFTW_MEASURE);
+    fft->reverse_plan = fftw_plan_dft_c2r_1d(nfft,
+            fft->complex_storage,
+            fft->real_storage,
+            FFTW_MEASURE);
+
+#endif
 
     /* print_vd(&fft->fft_work); */
-    
+
     feTRACE(15); 
     return fft;
 }
 void Fft_free(Fft* fft) {
     fsTRACE(15);
     dbgassert(fft,NULLPTRDEREF);
+#ifdef LASP_FFT_BACKEND_FFTPACK
     vd_free(&fft->fft_work);
+#elif defined LASP_FFT_BACKEND_FFTW
+    fftw_free(fft->complex_storage);
+    fftw_free(fft->real_storage);
+    fftw_destroy_plan(fft->forward_plan);
+    fftw_destroy_plan(fft->reverse_plan);
+#endif
     a_free(fft);
     feTRACE(15);
 }
+
 us Fft_nfft(const Fft* fft) {return fft->nfft;}
+
 void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) {
     fsTRACE(15);
     dbgassert(fft && freqdata && result,NULLPTRDEREF);
     const us nfft = fft->nfft;
     dbgassert(result->n_rows == nfft,
-              "Invalid size for time data rows."
-              " Should be equal to nfft");
+            "Invalid size for time data rows."
+            " Should be equal to nfft");
 
     dbgassert(freqdata->n_rows == (nfft/2+1),"Invalid number of rows in"
-              " result array");
+            " result array");
 
 
-    d* freqdata_ptr = (d*) getvcval(freqdata,0);
     d* result_ptr = getvdval(result,0);
 
+#ifdef LASP_FFT_BACKEND_FFTPACK
+    d* freqdata_ptr = (d*) getvcval(freqdata,0);
     /* Copy freqdata, to fft_result. */
     d_copy(&result_ptr[1],&freqdata_ptr[2],nfft-1,1,1);
     result_ptr[0] = freqdata_ptr[0];
 
     /* Perform inplace backward transform */
     npy_rfftb(nfft,
-              result_ptr,
-              getvdval(&fft->fft_work,0));
+            result_ptr,
+            getvdval(&fft->fft_work,0));
+
+
+#elif defined LASP_FFT_BACKEND_FFTW
+    c* freqdata_ptr = (c*) getvcval(freqdata,0);
+
+    c_copy(fft->complex_storage, freqdata_ptr,nfft/2+1);
+
+    fftw_execute(fft->reverse_plan);
+    
+    d_copy(result_ptr, fft->real_storage, nfft, 1, 1);
+    
+#endif
+    check_overflow_vx(*result);
 
     /* Scale by dividing by nfft. Checked with numpy implementation
      * that this indeed needs to be done for FFTpack. */
     d_scale(result_ptr,1/((d) nfft),nfft);
-    
-    check_overflow_vx(*result);
     feTRACE(15);
 }
 void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata) {
     fsTRACE(15);
-    
+
     dbgassert(fft && timedata && freqdata,NULLPTRDEREF);
 
     const us nchannels = timedata->n_cols;
     dbgassert(timedata->n_cols == freqdata->n_cols,
-              "Number of columns in timedata and result"
-              " should be equal.");
+            "Number of columns in timedata and result"
+            " should be equal.");
 
     for(us col=0;col<nchannels;col++) {
 
@@ -108,7 +151,7 @@ void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata) {
 }
 
 void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
-    
+
     fsTRACE(15);
     dbgassert(fft && timedata && result,NULLPTRDEREF);
 
@@ -116,12 +159,14 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
     assert_vx(timedata);
     assert_vx(result);
     dbgassert(timedata->n_rows == nfft,
-              "Invalid size for time data rows."
-              " Should be equal to nfft");
+            "Invalid size for time data rows."
+            " Should be equal to nfft");
 
     dbgassert(result->n_rows == (nfft/2+1),"Invalid number of rows in"
-              " result array");
+            " result array");
 
+
+#ifdef LASP_FFT_BACKEND_FFTPACK
     d* result_ptr = (d*) getvcval(result,0);
 
     /* Fftpack stores the data a bit strange, the resulting array
@@ -136,35 +181,45 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
 
     /* Perform fft */
     npy_rfftf(nfft,&result_ptr[1],
-              getvdval(&fft->fft_work,0));
+            getvdval(&fft->fft_work,0));
 
     /* Set real part of DC component to first index of the rfft
      * routine */
     result_ptr[0] = result_ptr[1];
 
     result_ptr[1] = 0;        /* Set imaginary part of DC component
-                                 * to zero */
+                               * to zero */
 
     /* For an even fft, the imaginary part of the Nyquist frequency
      * bin equals zero.*/
     if(likely(nfft%2 == 0)) {
         result_ptr[nfft+1] = 0;
     }
-
-    check_overflow_vx(*result);
     check_overflow_vx(fft->fft_work);
+#elif defined LASP_FFT_BACKEND_FFTW
+
+    d* timedata_ptr = getvdval(timedata,0);
+    c* result_ptr = getvcval(result,0);
+    d_copy(fft->real_storage,timedata_ptr, nfft, 1, 1);
+
+    fftw_execute(fft->forward_plan);
+    
+    c_copy(result_ptr, fft->complex_storage, nfft/2+1);
+    
+#endif
+    check_overflow_vx(*result);
     feTRACE(15);
 
 }
 void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result) {
     fsTRACE(15);
-    
+
     dbgassert(fft && timedata && result,NULLPTRDEREF);
 
     const us nchannels = timedata->n_cols;
     dbgassert(timedata->n_cols == result->n_cols,
-              "Number of columns in timedata and result"
-              " should be equal.");
+            "Number of columns in timedata and result"
+            " should be equal.");
 
     for(us col=0;col<nchannels;col++) {
 

lasp/c/lasp_firfilterbank.c

@@ -3,17 +3,17 @@
 // Author: J.A. de Jong -ASCEE
 // 
 // Description:
-// FilterBank implementation.
+// Firfilterbank implementation.
 //////////////////////////////////////////////////////////////////////
 #define TRACERPLUS (-5)
-#include "lasp_filterbank.h"
+#include "lasp_firfilterbank.h"
 #include "lasp_fft.h"
 #include "lasp_dfifo.h"
 #include "lasp_tracer.h"
 #include "lasp_alg.h"
 #define FIFO_SIZE_MULT 2
 
-typedef struct FilterBank_s {
+typedef struct Firfilterbank_s {
     us nfft;
 
     us P_m_1;                   /**< Filter length minus one */
@@ -25,9 +25,9 @@ typedef struct FilterBank_s {
 
     Fft* fft;                   /* Handle to internal FFT-function */
 
-} FilterBank;
+} Firfilterbank;
 
-FilterBank* FilterBank_create(const dmat* h,
+Firfilterbank* Firfilterbank_create(const dmat* h,
                               const us nfft) {
 
     fsTRACE(15);
@@ -46,7 +46,7 @@ FilterBank* FilterBank_create(const dmat* h,
         return NULL;
     }
 
-    FilterBank* fb = a_malloc(sizeof(FilterBank));
+    Firfilterbank* fb = a_malloc(sizeof(Firfilterbank));
 
     fb->nfft = nfft;
     fb->P_m_1 = P-1;
@@ -77,7 +77,7 @@ FilterBank* FilterBank_create(const dmat* h,
     feTRACE(15);
     return fb;
 }
-void FilterBank_free(FilterBank* fb) {
+void Firfilterbank_free(Firfilterbank* fb) {
     fsTRACE(15);
     dbgassert(fb,NULLPTRDEREF);
     cmat_free(&fb->filters);
@@ -87,7 +87,7 @@ void FilterBank_free(FilterBank* fb) {
     a_free(fb);
     feTRACE(15);
 }
-dmat FilterBank_filter(FilterBank* fb,
+dmat Firfilterbank_filter(Firfilterbank* fb,
                        const vd* x) {
 
     fsTRACE(15);

lasp/c/lasp_firfilterbank.h

@@ -1,8 +1,8 @@
-// lasp_filterbank.h
+// lasp_firfilterbank.h
 //
 // Author: J.A. de Jong - ASCEE
 //
-// Description: Implemententation of a discrete filterbank using fast
+// Description: Implemententation of a discrete FIR filterbank using fast
 // convolution and the overlap-save (overlap-scrap method). Multiple
 // filters can be applied to the same input data (*filterbank*).
 // Implementation is computationally efficient, as the forward FFT is
@@ -10,14 +10,14 @@
 // each filter in the filterbank.
 //////////////////////////////////////////////////////////////////////
 #pragma once
-#ifndef LASP_FILTERBANK_H
-#define LASP_FILTERBANK_H
+#ifndef LASP_FIRFILTERBANK_H
+#define LASP_FIRFILTERBANK_H
 #include "lasp_types.h"
 #include "lasp_mat.h"
-typedef struct FilterBank_s FilterBank;
+typedef struct Firfilterbank_s Firfilterbank;
 
 /** 
- * Initializes a fast convolution filter bank and returns a FilterBank
+ * Initializes a fast convolution filter bank and returns a Firfilterbank
  * handle. The nfft will be chosen to be at least four times the
  * length of the FIR filters.
  *
@@ -30,9 +30,9 @@ typedef struct FilterBank_s FilterBank;
  * times the filter lengths. For the lowest possible latency, it is
  * better to set nfft at twice the filter length.
  *
- * @return FilterBank handle, NULL on error.
+ * @return Firfilterbank handle, NULL on error.
  */
-FilterBank* FilterBank_create(const dmat* h,const us nfft);
+Firfilterbank* Firfilterbank_create(const dmat* h,const us nfft);
 
 /** 
  * Filters x using h, returns y
@@ -44,7 +44,7 @@ FilterBank* FilterBank_create(const dmat* h,const us nfft);
  * filterbank. The number of output samples is equal to the number of
  * input samples in x.
  */
-dmat FilterBank_filter(FilterBank* fb,
+dmat Firfilterbank_filter(Firfilterbank* fb,
                        const vd* x);
 
 /** 
@@ -52,8 +52,8 @@ dmat FilterBank_filter(FilterBank* fb,
  *
  * @param f Filter handle
  */
-void FilterBank_free(FilterBank* f);
+void Firfilterbank_free(Firfilterbank* f);
 
 
-#endif // LASP_FILTERBANK_H
+#endif // LASP_FIRFILTERBANK_H
 //////////////////////////////////////////////////////////////////////

lasp/c/lasp_mat.h

@@ -10,6 +10,7 @@
 #define LASP_MAT_H
 #include "lasp_math_raw.h"
 #include "lasp_alloc.h"
+#include "lasp_assert.h"
 #include "lasp_tracer.h"
 #include "lasp_assert.h"
 
@@ -18,15 +19,16 @@ typedef struct {
     us n_rows;
     us n_cols;
     bool _foreign_data;
-    us stride;
+    us colstride;
     d* _data;
 } dmat;
+
 /// Dense matrix of complex floating point values
 typedef struct {
     us n_rows;
     us n_cols;
     bool _foreign_data;
-    us stride;
+    us colstride;
     c* _data;
 } cmat;
 
@@ -46,9 +48,9 @@ typedef cmat vc;
     (vec)->_data[index] = val;
 
 #define setmatval(mat,row,col,val)                              \
-    dbgassert((((us) row) <= mat->n_rows),OUTOFBOUNDSMATR);     \
-    dbgassert((((us) col) <= mat->n_cols),,OUTOFBOUNDSMATC);    \
-    (mat)->data[(col)*(mat)->stride+(row)] = val;
+    dbgassert((((us) row) <= (mat)->n_rows),OUTOFBOUNDSMATR);     \
+    dbgassert((((us) col) <= (mat)->n_cols),OUTOFBOUNDSMATC);    \
+    (mat)->_data[(col)*(mat)->colstride+(row)] = val;
 
 /**
  * Return pointer to a value from a vector
@@ -67,7 +69,7 @@ static inline d* getdmatval(const dmat* mat,us row,us col){
     dbgassert(mat,NULLPTRDEREF);
     dbgassert(row < mat->n_rows,OUTOFBOUNDSMATR);
     dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC);
-    return &mat->_data[col*mat->stride+row];
+    return &mat->_data[col*mat->colstride+row];
 }
 /**
  * Return a value from a matrix of complex floating points
@@ -80,7 +82,7 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){
     dbgassert(mat,NULLPTRDEREF);
     dbgassert(row < mat->n_rows,OUTOFBOUNDSMATR);
     dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC);
-    return &mat->_data[col*mat->stride+row];
+    return &mat->_data[col*mat->colstride+row];
 }
 
 static inline d* getvdval(const vd* vec,us row){
@@ -109,7 +111,7 @@ static inline c* getvcval(const vc* vec,us row){
 #define check_overflow_xmat(xmat)                               \
     TRACE(15,"Checking overflow " #xmat);                        \
     if(!(xmat)._foreign_data) {                                  \
-        dbgassert((xmat)._data[((xmat).n_cols-1)*(xmat).stride+(xmat).n_rows] \
+        dbgassert((xmat)._data[((xmat).n_cols-1)*(xmat).colstride+(xmat).n_rows] \
                   == OVERFLOW_MAGIC_NUMBER,                             \
                   "Buffer overflow detected on" #xmat );                \
     }                                                                   \
@@ -164,7 +166,10 @@ static inline void cmat_set(cmat* mat,const c value){
  */
 static inline dmat dmat_alloc(us n_rows,
                               us n_cols) {
-    dmat result = { n_rows, n_cols, false, n_rows, NULL};
+    dmat result = { n_rows, n_cols, 
+                   false,
+                   n_rows, // The column stride
+                   NULL};
 
     #ifdef LASP_DEBUG
     result._data = (d*) a_malloc((n_rows*n_cols+1)*sizeof(d));
@@ -233,7 +238,7 @@ static inline vc vc_alloc(us size) {
  * Creates a dmat from foreign data. Does not copy the data, but only
  * initializes the row pointers. Assumes column-major ordering for the
  * data. Please do not keep this one alive after the data has been
- * destroyed. Assumes the column stride equals to n_rows.
+ * destroyed. Assumes the column colstride equals to n_rows.
  *
  * @param n_rows Number of rows
  * @param n_cols Number of columns
@@ -246,12 +251,12 @@ static inline dmat dmat_foreign(dmat* other) {
     dmat result = {other->n_rows,
                    other->n_cols,
                    true,
-                   other->stride,
+                   other->colstride,
                    other->_data};
     return result;
 }
 /**
- * Create a dmat from foreign data. Assumes the stride of the data is
+ * Create a dmat from foreign data. Assumes the colstride of the data is
  * n_rows.
  *
  * @param n_rows Number of rows
@@ -274,7 +279,7 @@ static inline dmat dmat_foreign_data(us n_rows,
     return result;
 }
 /**
- * Create a cmat from foreign data. Assumes the stride of the data is
+ * Create a cmat from foreign data. Assumes the colstride of the data is
  * n_rows.
  *
  * @param n_rows Number of rows
@@ -301,7 +306,7 @@ static inline cmat cmat_foreign_data(us n_rows,
  * Creates a cmat from foreign data. Does not copy the data, but only
  * initializes the row pointers. Assumes column-major ordering for the
  * data. Please do not keep this one alive after the data has been
- * destroyed. Assumes the column stride equals to n_rows.
+ * destroyed. Assumes the column colstride equals to n_rows.
  *
  * @param n_rows
  * @param n_cols
@@ -314,7 +319,7 @@ static inline cmat cmat_foreign(cmat* other) {
     cmat result = {other->n_rows,
                    other->n_cols,
                    true,
-                   other->stride,
+                   other->colstride,
                    other->_data};
     return result;
 }
@@ -391,7 +396,7 @@ static inline dmat dmat_submat(const dmat* parent,
 
     dmat result = { n_rows,n_cols,
                     true,           // Foreign data = true
-                    parent->n_rows, // This is the stride to get to
+                    parent->n_rows, // This is the colstride to get to
                     // the next column.
                     getdmatval(parent,startrow,startcol)};
 
@@ -421,7 +426,7 @@ static inline cmat cmat_submat(cmat* parent,
 
     cmat result = { n_rows,n_cols,
                     true,           // Foreign data = true
-                    parent->n_rows, // This is the stride to get to
+                    parent->n_rows, // This is the colstride to get to
                     // the next column.
                     getcmatval(parent,startrow,startcol)};
 
@@ -496,7 +501,7 @@ static inline void vc_copy(vc* to,const vc* from) {
  * @return vector with reference to column
  */
 static inline vd dmat_column(dmat* x,us col) {
-    vd res = {x->n_rows,1,true,x->stride,getdmatval(x,0,col)};
+    vd res = {x->n_rows,1,true,x->colstride,getdmatval(x,0,col)};
     return res;
 }
 
@@ -509,7 +514,7 @@ static inline vd dmat_column(dmat* x,us col) {
  * @return vector with reference to column
  */
 static inline vc cmat_column(cmat* x,us col) {
-    vc res = {x->n_rows,1,true,x->stride,getcmatval(x,0,col)};
+    vc res = {x->n_rows,1,true,x->colstride,getcmatval(x,0,col)};
     return res;
 }
 

lasp/c/lasp_math_raw.h

@@ -173,6 +173,8 @@ static inline d d_dot(const d a[],const d b[],const us size){
  * @param to : Array to write to
  * @param from : Array to read from
  * @param size : Size of arrays
+ * @param to_inc : Mostly equal to 1, the stride of the array to copy to
+ * @param from_inc : Mostly equal to 1, the stride of the array to copy from
  */
 static inline void d_copy(d to[],
                           const d from[],

lasp/c/lasp_python.h

@@ -20,10 +20,11 @@
 /** 
  * Create a numpy array from an existing dmat.
  *
- * @param mat 
- * @param transfer_ownership 
+ * @param mat dmat struccture containing array data and metadata.
+ * @param transfer_ownership If set to true, Numpy array will be responsible
+ * for freeing the data.
  *
- * @return 
+ * @return Numpy array
  */
 static inline PyObject* dmat_to_ndarray(dmat* mat,bool transfer_ownership) {
     fsTRACE(15);
@@ -50,15 +51,14 @@ static inline PyObject* dmat_to_ndarray(dmat* mat,bool transfer_ownership) {
     }
 
     // Transpose the array
-    PyObject* arr = PyArray_Transpose(arr_t,NULL);
+    PyObject* arr = PyArray_Transpose((PyArrayObject*) arr_t,NULL);
     if(!arr) {
         WARN("Array transpose failure");
         feTRACE(15);
         return NULL;
     }
     Py_DECREF(arr_t);
-
-                                              
+ 
     feTRACE(15);
     return arr;
 }

lasp/c/lasp_siggen.c

@@ -0,0 +1,182 @@
+// 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"
+
+#define PRIVATE_SIZE 64
+
+typedef enum {
+    SINEWAVE = 0,
+    WHITENOISE,
+    SWEEP,
+
+} SignalType;
+
+typedef struct Siggen {
+    SignalType signaltype;
+    d fs; // Sampling frequency [Hz]
+    d level_amp;
+    void* private;
+
+    char private_data[PRIVATE_SIZE];
+
+} Siggen;
+
+typedef struct { 
+    d curtime;
+    d omg;
+} SinewaveSpecific;
+
+typedef struct {
+    d V1, V2, S;
+    int phase;
+} WhitenoiseSpecific;
+
+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->private = NULL;
+    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* sinesiggen = Siggen_create(SINEWAVE, fs, level_dB);
+    sinesiggen->private = sinesiggen->private_data;
+    SinewaveSpecific* sp = (SinewaveSpecific*) sinesiggen->private;
+    sp->curtime = 0;
+    sp->omg = 2*number_pi*freq;
+
+    feTRACE(15);
+    return sinesiggen;
+}
+
+Siggen* Siggen_Whitenoise_create(const d fs, const d level_dB) {
+    fsTRACE(15);
+
+    Siggen* whitenoise = Siggen_create(WHITENOISE, fs, level_dB);
+    whitenoise->private = whitenoise->private_data;
+    dbgassert(sizeof(WhitenoiseSpecific) <= sizeof(whitenoise->private_data), "Allocated memory too small");
+    WhitenoiseSpecific* wn = whitenoise->private;
+    wn->phase = 0;
+    wn->V1 = 0;
+    wn->V2 = 0;
+    wn->S = 0;
+
+    feTRACE(15);
+    return whitenoise;
+}
+
+void Siggen_free(Siggen* siggen) {
+    fsTRACE(15);
+    assertvalidptr(siggen);
+
+    if(siggen->signaltype == SWEEP) {
+        /* Sweep specific stuff here */
+    }
+
+    a_free(siggen);
+    feTRACE(15);
+}
+
+static void Sinewave_genSignal(Siggen* siggen, SinewaveSpecific* sine, vd* samples) {
+    fsTRACE(15);
+    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(15);
+}
+
+static void Whitenoise_genSignal(Siggen* siggen, WhitenoiseSpecific* wn, vd* samples) {
+    fsTRACE(15);
+    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 * log(S) / S);
+        } else
+            X = V2 * sqrt(-2 * log(S) / S);
+
+        phase = 1 - phase;
+
+        setvecval(samples, i, siggen->level_amp*X);
+    }
+    wn->S = S;
+    wn->V1 = V1;
+    wn->V2 = V2;
+    wn->phase = phase;
+    feTRACE(15);
+
+}
+
+void Siggen_genSignal(Siggen* siggen,vd* samples) {
+
+    fsTRACE(15);
+    assertvalidptr(siggen);
+    assert_vx(samples);
+    d fs = siggen->fs;
+
+    switch(siggen->signaltype) {
+        case SINEWAVE:
+            Sinewave_genSignal(siggen, (SinewaveSpecific*) siggen->private, 
+                    samples);
+
+            break;
+        case WHITENOISE:
+            Whitenoise_genSignal(siggen, (WhitenoiseSpecific*) siggen->private,
+                    samples);
+            break;
+        case SWEEP:
+            break;
+        default:
+            dbgassert(false, "Not implementend signal type");
+
+    }
+
+    feTRACE(15);
+}
+
+
+
+//////////////////////////////////////////////////////////////////////

lasp/c/lasp_siggen.h

@@ -0,0 +1,58 @@
+// 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"
+
+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 white noise signal generator
+ * 
+ * @return Siggen* handle
+ */
+Siggen* Siggen_Whitenoise_create(const d fs, const d level_dB);
+
+/**
+ * Create a pink (1/f) noise signal generator
+ * 
+ * @param[in] fs: Sampling frequency [Hz]
+ * @return Siggen* handle
+ */
+Siggen* Siggen_Pinknoise_create(const us fs,const d level_dB);
+
+/* Siggen* Siggen_ForwardSweep_create(const d fs,; */
+/* Siggen* Siggen_(const d fs,; */
+
+/**
+ * Obtain a new piece of signal
+ * 
+ * @param[in] Siggen* Signal generator private data
+ * @param[out] samples Samples to fill. Vector should be pre-allocated
+ */
+void Siggen_genSignal(Siggen*,vd* samples);
+/**
+ * Free Siggen data
+ * 
+ * @param[in] Siggen* Signal generator private data
+ */
+void Siggen_free(Siggen*);
+
+#endif //LASP_SIGGEN_H
+//////////////////////////////////////////////////////////////////////

lasp/c/lasp_sosfilterbank.c

@@ -0,0 +1,150 @@
+#define TRACERPLUS 10
+#include "lasp_sosfilterbank.h"
+
+
+typedef struct Sosfilterbank {
+    /// The filter_coefs matrix contains filter coefficients for a SOS filter.
+    us filterbank_size;
+    us nsections;
+
+    /// The filter coefficients for each of the filters in the Filterbank
+    /// The *first* axis is the filter no, the second axis contains the
+    /// filter coefficients, in the order, b_0, b_1, b_2, a_0, a_1, a_2, which
+    /// corresponds to the transfer function 
+    ///                 b_0 + b_1 z^-1 + b_2 z^-2
+    ///         H[z] =  -------------------------
+    ///                 a_0 + a_1 z^-1 + a_2 z^-2 
+    dmat sos; /// sos[filter_no, coeff]
+
+    /// Storage for the current state of the output, first axis correspond to
+    /// the filter number axis, the second axis contains state coefficients 
+    dmat state;
+} Sosfilterbank;
+
+
+Sosfilterbank* Sosfilterbank_create(const us filterbank_size,
+        const us nsections) {
+    fsTRACE(15);
+    dbgassert(filterbank_size <= MAX_SOS_FILTER_BANK_SIZE, 
+            "Illegal filterbank size. Max size is " 
+            annestr(MAX_SOS_FILTER_BANK_SIZE));
+
+    Sosfilterbank* fb = (Sosfilterbank*) a_malloc(sizeof(Sosfilterbank));
+
+    fb->filterbank_size = filterbank_size;
+    dbgassert(nsections < MAX_SOS_FILTER_BANK_NSECTIONS,"Illegal number of sections");
+    fb->nsections = nsections;
+
+    /// Allocate filter coefficients matrix
+    fb->sos = dmat_alloc(filterbank_size, nsections*6);
+    fb->state = dmat_alloc(filterbank_size, nsections*2);
+    dmat_set(&(fb->state), 0);
+
+    /// Set all filter coefficients to unit impulse response
+    vd imp_response = vd_alloc(6*nsections);
+    vd_set(&imp_response,0);
+    for(us section = 0;section < nsections; section++) {
+        // Set b0 coefficient to 1
+        setvecval(&imp_response, 0 + 6*section, 1);
+        // Set a0 coefficient to 1
+        setvecval(&imp_response, 3 + 6*section, 1);
+    }
+
+    // Initialize all filters with a simple impulse response, single pass
+    for(us filter_no = 0; filter_no < filterbank_size; filter_no++) {
+        Sosfilterbank_setFilter(fb,filter_no,imp_response);
+    }
+    // Check if coefficients are properly initialized
+    // print_dmat(&(fb->sos));
+    vd_free(&imp_response);
+    feTRACE(15);
+    return fb;
+}
+
+void Sosfilterbank_setFilter(Sosfilterbank* fb,const us filter_no,
+        const vd filter_coefs) {
+    fsTRACE(15);
+    assertvalidptr(fb);
+    assert_vx(&filter_coefs);
+    dbgassert(filter_no < fb->filterbank_size, "Illegal filter number");
+    dbgassert(filter_coefs.n_rows == fb->nsections * 6, 
+            "Illegal filter coefficient length");
+
+    dmat *sos = &fb->sos;
+    dmat *state = &fb->state;
+    us nsections = fb->nsections;
+
+    for(us index=0;index<nsections*6;index++){
+        // Copy contents to position in sos matrix
+        *getdmatval(sos,filter_no,index) = *getvdval(&filter_coefs,index);
+    }
+
+    feTRACE(15);
+}
+
+void Sosfilterbank_free(Sosfilterbank* fb) {
+    fsTRACE(15);
+    assertvalidptr(fb);
+
+    dmat_free(&(fb->sos));
+    dmat_free(&(fb->state));
+
+    a_free(fb);
+    feTRACE(15);
+
+}
+
+dmat Sosfilterbank_filter(Sosfilterbank* fb,const vd* xs) {
+
+    fsTRACE(15);
+    assertvalidptr(fb);
+    assert_vx(xs);
+    dmat state = fb->state;
+    dmat sos = fb->sos;
+
+    us nsections = fb->nsections;
+    us filterbank_size = fb->filterbank_size;
+    us nsamples = xs->n_rows;
+
+    dmat ys = dmat_alloc(nsamples, filterbank_size);
+    /// Copy input signal to output array
+    for(us filter=0;filter<filterbank_size;filter++) {
+        d_copy(getdmatval(&ys,0,filter),getvdval(xs,0),nsamples,1,1);
+    }
+
+    /// Implementation is based on Proakis & Manolakis - Digital Signal
+    /// Processing, Fourth Edition, p. 550 
+    for(us section=0;section<nsections;section++) {
+        /// Obtain state information for current section, and all filters
+
+        for(us filter=0;filter<filterbank_size;filter++) {
+            d w1 = *getdmatval(&state,filter,section*2); 
+            d w2 = *getdmatval(&state,filter,section*2+1); 
+
+            d b0 = *getdmatval(&sos,filter,section*6+0);
+            d b1 = *getdmatval(&sos,filter,section*6+1);
+            d b2 = *getdmatval(&sos,filter,section*6+2);
+            d a0 = *getdmatval(&sos,filter,section*6+3);
+            d a1 = *getdmatval(&sos,filter,section*6+4);
+            d a2 = *getdmatval(&sos,filter,section*6+5);
+
+            d* y = getdmatval(&ys, 0, filter);
+
+            for(us sample=0;sample<nsamples;sample++){
+                d w0 = *y - a1*w1 - a2*w2;
+                d yn = b0*w0 + b1*w1 + b2*w2;
+                w2 = w1;
+                w1 = w0;
+                *y++ = yn;
+            }
+            *getdmatval(&state,filter,section*2)   = w1;
+            *getdmatval(&state,filter,section*2+1) = w2; 
+
+        }
+    }
+
+
+
+    feTRACE(15);
+    return ys;
+}

lasp/c/lasp_sosfilterbank.h

@@ -0,0 +1,62 @@
+// lasp_sosfilterbank.h
+//
+// Author: J.A. de Jong - ASCEE
+//
+// Description: Implemententation of a discrete filterbank using cascaded
+// second order sections (sos), also called BiQuads.
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef LASP_FILTERBANK_H
+#define LASP_FILTERBANK_H
+#include "lasp_types.h"
+#include "lasp_mat.h"
+
+#define MAX_SOS_FILTER_BANK_SIZE 40
+#define MAX_SOS_FILTER_BANK_NSECTIONS 6
+
+typedef struct Sosfilterbank Sosfilterbank;
+
+/** 
+ * Initializes a Sosfilterbank. Sets all coefficients in such a way that the
+ * filter effectively does nothing (unit impulse response).
+ */
+Sosfilterbank* Sosfilterbank_create(const us filterbank_size,
+                                    const us nsections);
+
+/**
+ * Initialize the filter coeficients in the filterbank
+ *
+ * @param fb: Filterbank handle
+ * @param filter_no: Filter number in the bank
+ * @param coefss: Array of filter coefficients. Should have a length of 
+ * nsections x 6, for each of the sections, it contains (b0, b1, b2, a0, 
+ * a1, a2), where a are the numerator coefficients and b are the denominator 
+ * coefficients. 
+ *
+*/
+void Sosfilterbank_setFilter(Sosfilterbank* fb,const us filter_no,
+                             const vd coefs);
+
+/** 
+ * Filters x using h, returns y
+ *
+ * @param x Input time sequence block. Should have at least one sample.
+
+ * @return Filtered output in an allocated array. The number of
+ * columns in this array equals the number of filters in the
+ * filterbank. The number of output samples is equal to the number of
+ * input samples in x.
+ */
+dmat Sosfilterbank_filter(Sosfilterbank* fb,
+                          const vd* x);
+
+/** 
+ * Cleans up an existing filter bank.
+ *
+ * @param f Filterbank handle
+ */
+void Sosfilterbank_free(Sosfilterbank* f);
+
+
+#endif // LASP_FILTERBANK_H
+//////////////////////////////////////////////////////////////////////

lasp/config.pxi.in

@@ -1,68 +1,10 @@
 import numpy as np
-cimport numpy as np
-from libcpp cimport bool
+cimport numpy as cnp
+
+# Do this, our segfaults will be your destination
+cnp.import_array()
 
 DEF LASP_FLOAT = "@LASP_FLOAT@"
+DEF LASP_DEBUG_CYTHON = "@LASP_DEBUG_CYTHON@"
 
-IF LASP_FLOAT == "double":
-    ctypedef  double d
-    ctypedef double complex c
-    NUMPY_FLOAT_TYPE = np.float64
-    NUMPY_COMPLEX_TYPE = np.complex128
-    CYTHON_NUMPY_FLOAT_t = np.NPY_FLOAT64
-    CYTHON_NUMPY_COMPLEX_t = np.NPY_COMPLEX128
-
-ELSE:
-    ctypedef  float d
-    ctypedef float complex c
-    NUMPY_FLOAT_TYPE = np.float32
-    NUMPY_COMPLEX_TYPE = np.complex64
-    CYTHON_NUMPY_FLOAT_t = np.NPY_FLOAT32
-    CYTHON_NUMPY_COMPLEX_t = np.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":
-    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(np.ndarray arr, int flags)
-
-
-cdef extern from "lasp_python.h":
-    object dmat_to_ndarray(dmat*,bint transfer_ownership)
+from libcpp cimport bool

lasp/device/CMakeLists.txt

@@ -1,4 +1,9 @@
-set_source_files_properties(lasp_daqdevice.c PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} ${CYTHON_EXTRA_C_FLAGS}")
-cython_add_module(lasp_daqdevice lasp_daqdevice.pyx)
+include_directories(/usr/include/rtaudio)
+set_source_files_properties(lasp_rtaudio.pyx PROPERTIES CYTHON_IS_CXX TRUE)
+set_source_files_properties(lasp_rtaudio.cxx PROPERTIES COMPILE_FLAGS
+    "${CMAKE_CXX_FLAGS} ${CYTHON_EXTRA_CXX_FLAGS}")
 
-target_link_libraries(lasp_daqdevice asound)
+cython_add_module(lasp_rtaudio lasp_rtaudio.pyx)
+
+# target_link_libraries(lasp_portaudio portaudio) */
+target_link_libraries(lasp_rtaudio rtaudio)

lasp/device/__init__.py

@@ -1,6 +1,6 @@
 #!/usr/bin/python3
-__all__ = ['DAQDevice', 'DAQConfiguration', 'configs', 'query_devices',
-           'roga_plugndaq', 'default_soundcard']
-from .lasp_daqdevice import DAQDevice, query_devices
-from .lasp_daqconfig import (DAQConfiguration, configs, roga_plugndaq,
-                             default_soundcard)
+__all__ = ['DAQConfiguration']
+from .lasp_daqconfig import DAQConfiguration, DAQInputChannel, DeviceInfo
+from .lasp_rtaudio import (RtAudio, 
+                           get_numpy_dtype_from_format_string,
+                           get_sampwidth_from_format_string)

lasp/device/lasp_daqconfig.py

@@ -9,10 +9,31 @@ Data Acquistiion (DAQ) device descriptors, and the DAQ devices themselves
 
 """
 from dataclasses import dataclass, field
-from .lasp_daqdevice import query_devices, DeviceInfo
-__all__ = ['DAQConfiguration', 'roga_plugndaq', 'umik',
-           'default_soundcard', 'configs',
-           'findDAQDevice']
+import numpy as np
+
+
+@dataclass
+class DeviceInfo:
+    index: int
+    probed: bool
+    name: str
+    outputchannels: int
+    inputchannels: int
+    duplexchannels: int
+    samplerates: list
+    sampleformats: list
+    prefsamplerate: int
+
+
+from ..lasp_common import lasp_shelve
+
+@dataclass
+class DAQInputChannel:
+    channel_enabled: bool
+    channel_name: str
+    sensitivity: float
+
+
 
 
 @dataclass
@@ -21,11 +42,11 @@ class DAQConfiguration:
     Initialize a device descriptor
 
     Args:
-        name: Name of the device to appear to the user
-        cardname: ALSA name identifier
-        cardlongnamematch: Long name according to ALSA
-        device_name: ASCII name with which to open the device when connected
-        en_format: index in the list of sample formats
+        input_device_name: ASCII name with which to open the device when connected
+        outut_device_name: ASCII name with which to open the device when connected
+
+        ==============================
+        en_format: index of the format in the list of sample formats
         en_input_rate: index of enabled input sampling frequency [Hz]
                             in the list of frequencies.
         en_input_channels: list of channel indices which are used to
@@ -39,105 +60,55 @@ class DAQConfiguration:
         en_output_rate: index in the list of possible output sampling
                              frequencies.
         en_output_channels: list of enabled output channels
+        ===============================
 
 
     """
+    duplex_mode: bool
 
-    name: str
-    cardname: str
-    cardlongnamematch: str
-    device_name: str
-    en_format: int
+    input_device_name: str
+    output_device_name: str
+
+    en_input_sample_format: str
+    en_output_sample_format: str
     en_input_rate: int
+    en_output_rate: int
+
     en_input_channels: list
-    input_sensitivity: list
-    input_gain_settings: list
-    en_input_gain_settings: list = field(default_factory=list)
-    en_output_rate: int = -1
-    en_output_channels: list = field(default_factory=list)
 
-    def match(self, device):
+    def getEnabledChannels(self):
+        en_channels = []
+        for chan in self.en_input_channels:
+            if chan.channel_enabled:
+                en_channels.append(chan)
+        return en_channels
+
+    def getSensitivities(self):
+        return np.array(
+                [float(channel.sensitivity) for channel in
+                    self.getEnabledChannels()])
+
+    @staticmethod
+    def loadConfigs():
         """
-        Returns True when a device specification dictionary matches to the
-        configuration.
-
-        Args:
-            device: dictionary specifying device settings
+        Returns a list of currently available configurations
         """
-        match = True
-        if self.cardlongnamematch is not None:
-            match &= self.cardlongnamematch in device.cardlongname
-        if self.cardname is not None:
-            match &= self.cardname == device.cardname
-        if self.device_name is not None:
-            match &= self.device_name == device.device_name
-        match &= self.en_format < len(device.available_formats)
-        match &= self.en_input_rate < len(device.available_input_rates)
-        match &= max(
-            self.en_input_channels) < device.max_input_channels
-        if len(self.en_output_channels) > 0:
-            match &= max(
-                self.en_output_channels) < device.max_output_channels
-        match &= self.en_output_rate < len(
-            device.available_output_rates)
+        with lasp_shelve() as sh:
+            return sh['daqconfigs'] if 'daqconfigs' in sh.keys() else {}
 
-        return match
+    def saveConfig(self, name):
+        with lasp_shelve() as sh:
+            if 'daqconfigs' in sh.keys():
+                cur_configs = sh['daqconfigs']
+            else:
+                cur_configs = {}
+            cur_configs[name] = self
+            sh['daqconfigs'] = cur_configs
 
+    @staticmethod
+    def deleteConfig(name):
+        with lasp_shelve() as sh:
+            cur_configs = sh['daqconfigs']
+            del cur_configs[name]
+            sh['daqconfigs'] = cur_configs
 
-roga_plugndaq = DAQConfiguration(name='Roga-instruments Plug.n.DAQ USB',
-                                 cardname='USB Audio CODEC',
-                                 cardlongnamematch='Burr-Brown from TI USB'
-                                 ' Audio CODEC',
-                                 device_name='iec958:CARD=CODEC,DEV=0',
-                                 en_format=0,
-                                 en_input_rate=2,
-                                 en_input_channels=[0],
-                                 input_sensitivity=[46.92e-3, 46.92e-3],
-                                 input_gain_settings=[-20, 0, 20],
-                                 en_input_gain_settings=[1, 1],
-                                 en_output_rate=1,
-                                 en_output_channels=[False, False]
-                                 )
-umik = DAQConfiguration(name='UMIK-1',
-                        cardname='Umik-1  Gain: 18dB',
-                        cardlongnamematch='miniDSP Umik-1  Gain: 18dB',
-                        device_name='iec958:CARD=U18dB,DEV=0',
-                        en_format=0,
-                        en_input_rate=0,
-                        en_input_channels=[0],
-                        input_sensitivity=[1., 1.],
-                        input_gain_settings=[0., 0.],
-                        en_input_gain_settings=[0, 0],
-                        en_output_rate=0,
-                        en_output_channels=[True, True]
-                        )
-
-default_soundcard = DAQConfiguration(name="Default device",
-                                     cardname=None,
-                                     cardlongnamematch=None,
-                                     device_name='default',
-                                     en_format=0,
-                                     en_input_rate=2,
-                                     en_input_channels=[0],
-                                     input_sensitivity=[1.0, 1.0],
-                                     input_gain_settings=[0],
-                                     en_input_gain_settings=[0, 0],
-                                     en_output_rate=1,
-                                     en_output_channels=[]
-                                     )
-configs = (roga_plugndaq, default_soundcard)
-
-
-def findDAQDevice(config: DAQConfiguration) -> DeviceInfo:
-    """
-    Search for a DaQ device for the given configuration.
-
-    Args:
-        config: configuration to search a device for
-    """
-    devices = query_devices()
-
-    for device in devices:
-        if config.match(device):
-            return device
-    return None

lasp/device/lasp_daqdevice.pyx

@@ -1,504 +0,0 @@
-include "config.pxi"
-from libc.stdlib cimport malloc, free
-from libc.stdio cimport printf, stderr, fprintf
-import sys
-import numpy as np
-cimport numpy as cnp
-
-__all__ = ['DAQDevice']
-
-from libc.errno cimport EPIPE, EBADFD, ESTRPIPE
-
-
-cdef extern from "alsa/asoundlib.h":
-    int snd_card_get_longname(int index,char** name)
-    int snd_card_get_name(int index,char** name)
-    int snd_card_next(int* rcard)
-
-    ctypedef struct snd_pcm_t
-    ctypedef struct snd_pcm_info_t
-    ctypedef struct snd_pcm_hw_params_t
-    ctypedef enum snd_pcm_stream_t:
-        SND_PCM_STREAM_PLAYBACK
-        SND_PCM_STREAM_CAPTURE
-    ctypedef enum snd_pcm_format_t:
-        SND_PCM_FORMAT_S16_LE
-        SND_PCM_FORMAT_S16_BE
-        SND_PCM_FORMAT_U16_LE
-        SND_PCM_FORMAT_U16_BE
-        SND_PCM_FORMAT_S24_LE
-        SND_PCM_FORMAT_S24_BE
-        SND_PCM_FORMAT_U24_LE
-        SND_PCM_FORMAT_U24_BE
-        SND_PCM_FORMAT_S32_LE
-        SND_PCM_FORMAT_S32_BE
-        SND_PCM_FORMAT_U32_LE
-        SND_PCM_FORMAT_U32_BE
-        SND_PCM_FORMAT_S24_3LE
-        SND_PCM_FORMAT_S24_3BE
-        SND_PCM_FORMAT_U24_3LE
-        SND_PCM_FORMAT_U24_3BE
-        SND_PCM_FORMAT_S16
-        SND_PCM_FORMAT_U16
-        SND_PCM_FORMAT_S24
-        SND_PCM_FORMAT_U24
-    const char* snd_pcm_format_name (snd_pcm_format_t format)
-    ctypedef enum snd_pcm_access_t:
-        SND_PCM_ACCESS_RW_INTERLEAVED
-    ctypedef unsigned long snd_pcm_uframes_t
-    int snd_pcm_open(snd_pcm_t** pcm,char* name, snd_pcm_stream_t type, int mode)
-    int snd_pcm_close(snd_pcm_t*)
-
-    int snd_pcm_hw_params_set_access(snd_pcm_t*, snd_pcm_hw_params_t*, snd_pcm_access_t)
-    void snd_pcm_hw_params_alloca(snd_pcm_hw_params_t**)
-    int snd_pcm_hw_params_any(snd_pcm_t*, snd_pcm_hw_params_t* params)
-    int snd_pcm_hw_params_current(snd_pcm_t*, snd_pcm_hw_params_t* params)
-    int snd_pcm_hw_params_set_rate_resample(snd_pcm_t*, snd_pcm_hw_params_t*, int)
-    int snd_pcm_hw_params_set_rate(snd_pcm_t*, snd_pcm_hw_params_t*,unsigned int val,int dir)
-    int snd_pcm_hw_params_set_format(snd_pcm_t*, snd_pcm_hw_params_t*, snd_pcm_format_t)
-    int snd_pcm_hw_params_set_channels(snd_pcm_t*, snd_pcm_hw_params_t*, unsigned val)
-    int snd_pcm_hw_params_set_period_size_near(snd_pcm_t*,snd_pcm_hw_params_t*,
-        snd_pcm_uframes_t*,int* dir)
-    int snd_pcm_hw_params_get_period_size(snd_pcm_hw_params_t*,
-        snd_pcm_uframes_t*,int* dir)
-    int snd_pcm_info(snd_pcm_t*, snd_pcm_info_t*)
-    void snd_pcm_info_alloca(snd_pcm_info_t**)
-    int snd_pcm_info_get_card(snd_pcm_info_t*)
-
-    int snd_pcm_drain(snd_pcm_t*)
-    int snd_pcm_readi(snd_pcm_t*,void* buf,snd_pcm_uframes_t nframes) nogil
-    int snd_pcm_hw_params(snd_pcm_t*,snd_pcm_hw_params_t*)
-    int snd_pcm_hw_params_test_rate(snd_pcm_t*, snd_pcm_hw_params_t*,
-                                    unsigned int val,int dir)
-    int snd_pcm_hw_params_test_format(snd_pcm_t*, snd_pcm_hw_params_t*,
-                                      snd_pcm_format_t)
-    int snd_pcm_hw_params_get_channels_max(snd_pcm_hw_params_t*,unsigned int*)
-    int snd_device_name_hint(int card, const char* iface, void*** hints)
-    char* snd_device_name_get_hint(void* hint, const char* id)
-    int snd_device_name_free_hint(void**)
-    char* snd_strerror(int rval)
-
-# Check for these sample rates
-check_rates = [8000, 44100, 48000, 96000, 19200]
-
-# First value in tuple: number of significant bits
-# Second value:         number of bits used in memory
-# Third value:          S for signed, U for unsigned, L for little endian,
-#                       and B for big endian.
-check_formats = {SND_PCM_FORMAT_S16_LE: (16,16,'SL'),
-                 SND_PCM_FORMAT_S16_BE: (16,16,'SB'),
-                 SND_PCM_FORMAT_U16_LE: (16,16,'UL'),
-                 SND_PCM_FORMAT_U16_BE: (16,16,'UB'),
-                 SND_PCM_FORMAT_S24_LE: (24,32,'SL'),
-                 SND_PCM_FORMAT_S24_BE: (24,32,'SB'),
-                 SND_PCM_FORMAT_U24_LE: (24,32,'UL'),
-                 SND_PCM_FORMAT_U24_BE: (24,32,'UB'),
-                 SND_PCM_FORMAT_S32_LE: (32,32,'SL'),
-                 SND_PCM_FORMAT_S32_BE: (32,32,'SB'),
-                 SND_PCM_FORMAT_U32_LE: (32,32,'UL'),
-                 SND_PCM_FORMAT_U32_BE: (32,32,'UB'),
-                 SND_PCM_FORMAT_S24_3LE: (24,24,'SL'),
-                 SND_PCM_FORMAT_S24_3BE: (24,24,'SB'),
-                 SND_PCM_FORMAT_U24_3LE: (24,24,'UL'),
-                 SND_PCM_FORMAT_U24_3BE: (24,24,'UB')}
-
-devices_opened_card = [False, False, False, False, False, False]
-
-cdef snd_pcm_t* open_device(char* name,
-                            snd_pcm_stream_t streamtype):
-    """
-    Helper function to properly open the first device of a card
-    """
-    cdef snd_pcm_t* pcm
-    # if name in devices_opened:
-        # raise RuntimeError('Device %s is already opened.' %name)
-    cdef int rval = snd_pcm_open(&pcm, name, streamtype, 0)
-    if rval == 0:
-        return pcm
-    else:
-        return NULL
-
-cdef int close_device(snd_pcm_t* dev):
-    rval = snd_pcm_close(dev)
-    if rval:
-        print('Error closing device')
-    return rval
-
-class DeviceInfo:
-    """
-    Will later be replaced by a dataclass. Storage container for a lot of
-    device parameters.
-    """
-
-    def __repr__(self):
-        rep = f"""Device name: {self.device_name}
-Card name: {self.cardname}
-Available sample formats: {self.available_formats}
-Max input channels: {self.max_input_channels}
-
-        """
-        return rep
-
-def getDeviceInfo(char* device_name):
-    """
-    Open the PCM device for both capture and playback, extract the number of
-    channels, the samplerates and encoding
-
-    Args:
-        cardindex: Card number of device, numbered by ALSA
-
-    Returns:
-
-    """
-    cdef:
-        snd_pcm_t* pcm
-        snd_pcm_hw_params_t* hwparams
-        snd_pcm_info_t* info
-        int rval, cardindex
-        unsigned max_input_channels=0, max_output_channels=0
-        char* c_cardname
-        char *c_cardlongname
-
-    deviceinfo = DeviceInfo()
-    deviceinfo.device_name = device_name.decode('ASCII')
-    pcm = open_device(device_name, SND_PCM_STREAM_CAPTURE)
-    if not pcm:
-        raise RuntimeError('Unable to open device')
-
-    snd_pcm_info_alloca(&info)
-    rval = snd_pcm_info(pcm, info)
-    if rval:
-        snd_pcm_close(pcm)
-        raise RuntimeError('Unable to obtain device info')
-
-    cardindex = snd_pcm_info_get_card(info)
-    cardname = ''
-    cardlongname = ''
-    if cardindex >= 0:
-        snd_card_get_name(cardindex, &c_cardname)
-        if c_cardname:
-            cardname = c_cardname.decode('ASCII')
-            printf('name: %s\n', c_cardname)
-            free(c_cardname)
-
-        rval = snd_card_get_longname(cardindex, &c_cardlongname)
-        if c_cardlongname:
-            printf('longname: %s\n', c_cardlongname)
-            cardlongname = c_cardlongname.decode('ASCII')
-            free(c_cardlongname)
-    deviceinfo.cardname = cardname
-    deviceinfo.cardlongname = cardlongname
-
-    # Check hardware parameters
-    snd_pcm_hw_params_alloca(&hwparams)
-
-    # Nothing said about the return value of this function in the documentation
-    snd_pcm_hw_params_any(pcm, hwparams)
-
-    # Check available sample formats
-    available_formats = []
-    for format in check_formats.keys():
-        rval = snd_pcm_hw_params_test_format(pcm, hwparams, format)
-        if rval == 0:
-            available_formats.append(check_formats[format])
-    deviceinfo.available_formats = available_formats
-    # # Restrict a configuration space to contain only real hardware rates.
-    # rval = snd_pcm_hw_params_set_rate_resample(pcm, hwparams, 0)
-    # if rval !=0:
-    #     fprintf(stderr, 'Unable disable resampling rates')
-
-    # Check available input sample rates
-    available_input_rates = []
-    for rate in check_rates:
-        rval = snd_pcm_hw_params_test_rate(pcm, hwparams, rate, 0)
-        if rval == 0:
-            available_input_rates.append(rate)
-    deviceinfo.available_input_rates = available_input_rates
-
-    rval = snd_pcm_hw_params_get_channels_max(hwparams, &max_input_channels)
-    if rval != 0:
-        fprintf(stderr, "Could not obtain max input channels\n")
-    deviceinfo.max_input_channels = max_input_channels
-
-    # Close device
-    rval = snd_pcm_close(pcm)
-    if rval:
-        fprintf(stderr, 'Unable to close pcm device.\n')
-
-    deviceinfo.available_output_rates = []
-    deviceinfo.max_output_channels = 0
-
-    # ###############################################################
-    # Open device for output
-    pcm = open_device(device_name, SND_PCM_STREAM_CAPTURE)
-    if pcm == NULL:
-        # We are unable to open the device for playback, but we were able to
-        # open in for capture. So this is a valid device.
-        return deviceinfo
-
-    snd_pcm_hw_params_any(pcm, hwparams)
-    # Restrict a configuration space to contain only real hardware rates.
-    # rval = snd_pcm_hw_params_set_rate_resample(pcm, hwparams, 0)
-    # if rval != 0:
-    #     fprintf(stderr, 'Unable disable resampling rates')
-
-    # Check available input sample rates
-    available_output_rates = []
-    for rate in check_rates:
-        rval = snd_pcm_hw_params_test_rate(pcm, hwparams, rate, 0)
-        if rval == 0:
-            available_output_rates.append(rate)
-    deviceinfo.available_output_rates = available_output_rates
-    rval = snd_pcm_hw_params_get_channels_max(hwparams, &max_output_channels)
-    if rval != 0:
-        fprintf(stderr, "Could not obtain max output channels")
-    deviceinfo.max_output_channels = max_output_channels
-
-    # Close device
-    rval = close_device(pcm)
-    if rval:
-        fprintf(stderr, 'Unable to close pcm device %s.', device_name)
-
-    return deviceinfo
-
-
-def query_devices():
-    """
-    Returns a list of available DAQ devices, where each device is represented
-    by a dictionary containing parameters of the device
-    """
-
-    devices = []
-
-    cdef:
-        # Start cardindex at -1, such that the first one is picked by
-        # snd_card_next()
-        int cardindex = -1, rval=0, i=0
-        void** namehints_opaque
-        char** namehints
-        char* c_device_name
-
-    rval = snd_device_name_hint(-1, "pcm", &namehints_opaque)
-    if rval:
-        raise RuntimeError('Could not obtain name hints for card %i.'
-        %cardindex)
-
-    namehints = <char**> namehints_opaque
-    while namehints[i] != NULL:
-        # printf('namehint[i]: %s\n', namehints[i])
-        c_device_name = snd_device_name_get_hint(namehints[i], "NAME")
-        c_descr = snd_device_name_get_hint(namehints[i], "DESC")
-        if c_device_name:
-            device_name = c_device_name.decode('ASCII')
-            if c_descr:
-                device_desc = c_descr.decode('ASCII')
-                free(c_descr)
-            else:
-                device_desc = ''
-            try:
-                device = getDeviceInfo(c_device_name)
-                printf('device name: %s\n', c_device_name)
-                devices.append(device)
-            except RuntimeError:
-                pass
-            free(c_device_name)
-
-        i+=1
-
-    snd_device_name_free_hint(namehints_opaque)
-
-    return devices
-
-
-
-cdef class DAQDevice:
-    cdef:
-        snd_pcm_t* pcm
-        int device_index
-        object device, config
-        public snd_pcm_uframes_t blocksize
-        object callback
-
-
-    def __cinit(self):
-        self.pcm = NULL
-
-    def __init__(self, config, blocksize=2048):
-        """
-        Initialize the DAQ device
-
-        Args:
-            config: DAQConfiguration instance
-            blocksize: Number of frames in a single acquisition block
-            callback: callback used to send data frames to
-        """
-
-        self.config = config
-        devices = query_devices()
-
-        self.device = None
-        for device in devices:
-            if self.config.match(device):
-                # To open the underlying PCM device
-                device_name = device.device_name.encode('ASCII')
-                self.device = device
-
-        if self.device is None:
-            raise RuntimeError(f'Device {self.config.name} is not available')
-        # if devices_opened[device_index]:
-            # raise RuntimeError(f'Device {self.config.name} is already opened')
-        # print('device_name opened:', device_name)
-        self.pcm = open_device(device_name, SND_PCM_STREAM_CAPTURE)
-
-        # Device is opened. We are going to configure
-        cdef:
-           snd_pcm_hw_params_t* params
-           int rval
-           snd_pcm_format_t format_code
-           snd_pcm_uframes_t period_size
-        snd_pcm_hw_params_alloca(&params);
-
-        # Fill it in with default values.
-        snd_pcm_hw_params_any(self.pcm, params);
-
-        # Set access interleaved
-        rval = snd_pcm_hw_params_set_access(self.pcm, params,
-         SND_PCM_ACCESS_RW_INTERLEAVED)
-        if rval != 0:
-            snd_pcm_close(self.pcm)
-            raise RuntimeError('Could not set access mode to interleaved')
-
-        # Set sampling frequency
-        cdef unsigned int rate
-        rate = device.available_input_rates[config.en_input_rate]
-        # printf('Set sample rate: %i\n', rate)
-        rval = snd_pcm_hw_params_set_rate(self.pcm,params, rate, 0)
-        if rval != 0:
-            snd_pcm_close(self.pcm)
-            raise RuntimeError('Could not set input sampling frequency')
-
-        # Set number of channels
-        channels_max = max(self.channels_en)+1
-        # print('channels_max:', channels_max)
-        if channels_max > self.device.max_input_channels:
-            snd_pcm_close(self.pcm)
-            raise ValueError('Highest required channel is larger than available'
-            ' channels.')
-        rval = snd_pcm_hw_params_set_channels(self.pcm, params, channels_max)
-        if rval != 0:
-            snd_pcm_close(self.pcm)
-            raise RuntimeError('Could not set input channels, highest required'
-            ' input channel: %i.' %channels_max)
-
-        # Find the format description
-        format_descr = self.device.available_formats[config.en_format]
-        # Obtain key from value of dictionary
-
-        format_code = list(check_formats.keys())[list(check_formats.values()).index(format_descr)]
-        # printf('Format code: %s\n', snd_pcm_format_name(format_code))
-
-        # print(format)
-        rval = snd_pcm_hw_params_set_format(self.pcm, params, format_code)
-        if rval != 0:
-            fprintf(stderr, "Could not set format: %s.", snd_strerror(rval))
-
-        # Set period size
-        cdef int dir = 0
-        bytedepth = format_descr[1]//8
-        # print('byte depth:', bytedepth)
-        period_size = blocksize
-        rval = snd_pcm_hw_params_set_period_size_near(self.pcm,
-                                                      params,
-                                                      &period_size, &dir)
-
-        if rval != 0:
-            snd_pcm_close(self.pcm)
-            raise RuntimeError("Could not set period size: %s."
-            %snd_strerror(rval))
-
-        # Write the parameters to the driver
-        rc = snd_pcm_hw_params(self.pcm, params);
-        if (rc < 0):
-            snd_pcm_close(self.pcm)
-            raise RuntimeError('Could not set hw parameters: %s' %snd_strerror(rc))
-
-        # Check the block size again, and store it
-        snd_pcm_hw_params_get_period_size(params, &self.blocksize,
-                                          &dir)
-        # print('Period size:', self.blocksize)
-
-    cdef object _getEmptyBuffer(self):
-        """
-        Return right size empty buffer
-        """
-        format_descr =  self.device.available_formats[self.config.en_format]
-        LB = format_descr[2][1]
-        assert LB == 'L', 'Only little-endian data format supported'
-        if format_descr[1] == 16:
-            dtype = np.int16
-        elif format_descr[1] == 32:
-            dtype = np.int32
-
-        # interleaved data, order = C
-        return np.zeros((self.blocksize,
-                         max(self.channels_en)+1),
-                          dtype=dtype, order='C')
-
-    def read(self):
-        cdef int rval = 0
-        buf = self._getEmptyBuffer()
-        # buf2 = self._getEmptyBuffer()
-        cdef cnp.int16_t[:, ::1] bufv = buf
-        with nogil:
-            rval = snd_pcm_readi(self.pcm,<void*> &bufv[0, 0], self.blocksize)
-        # rval = 2048
-        if rval > 0:
-            # print('Samples obtained:' , rval)
-            return buf[:rval, self.channels_en]
-
-            # return buf
-        elif rval == -EPIPE:
-            raise RuntimeError('Error: buffer overrun: %s',
-            snd_strerror(rval))
-        elif rval == -EBADFD:
-            raise RuntimeError('Error: could not read from DAQ Device: %s',
-            snd_strerror(rval))
-        elif rval == -ESTRPIPE:
-            raise RuntimeError('Error: could not read from DAQ Device: %s',
-            snd_strerror(rval))
-
-
-    def __dealloc__(self):
-        # printf("dealloc\n")
-        cdef int rval
-        if self.pcm:
-            # print('Closing pcm')
-            # snd_pcm_drain(self.pcm)
-            rval = snd_pcm_close(self.pcm)
-            # devices_opened[self.device_index] = False
-            if rval != 0:
-                fprintf(stderr, 'Unable to properly close device: %s\n',
-                snd_strerror(rval))
-            self.pcm = NULL
-
-    @property
-    def nchannels_all(self):
-        return self.device.max_input_channels
-
-    @property
-    def channels_en(self):
-        return self.config.en_input_channels
-
-    @property
-    def input_rate(self):
-        return self.device.available_input_rates[self.config.en_input_rate]
-
-    @property
-    def channels(self):
-        return [self.config.en_input_channels]
-
-    cpdef bint isOpened(self):
-        if self.pcm:
-            return True
-        else:
-            return False

lasp/device/lasp_rtaudio.pyx

@@ -0,0 +1,383 @@
+import sys
+include "config.pxi"
+cimport cython
+from .lasp_daqconfig import DeviceInfo
+from libcpp.string cimport string
+from libcpp.vector cimport vector
+from libc.string cimport memcpy
+
+cdef extern from "RtAudio.h" nogil:
+    ctypedef unsigned long RtAudioStreamStatus
+    RtAudioStreamStatus RTAUDIO_INPUT_OVERFLOW
+    RtAudioStreamStatus RTAUDIO_OUTPUT_UNDERFLOW
+
+    cdef cppclass RtAudioError:
+        ctypedef enum Type:
+            WARNING
+            DEBUG_WARNING
+            UNSPECIFIED
+            NO_DEVICES_FOUND
+            INVALID_DEVICE
+            MEMORY_ERROR
+            INVALID_PARAMETER
+            INVALID_USE
+            DRIVER_ERROR
+            SYSTEM_ERROR
+            THREAD_ERROR
+
+    ctypedef unsigned long RtAudioStreamFlags
+    RtAudioStreamFlags RTAUDIO_NONINTERLEAVED
+    RtAudioStreamFlags RTAUDIO_MINIMIZE_LATENCY
+    RtAudioStreamFlags RTAUDIO_HOG_DEVICE
+    RtAudioStreamFlags RTAUDIO_ALSA_USE_DEFAULT
+    RtAudioStreamFlags RTAUDIO_JACK_DONT_CONNECT
+
+    ctypedef unsigned long RtAudioFormat
+    RtAudioFormat RTAUDIO_SINT8
+    RtAudioFormat RTAUDIO_SINT16
+    RtAudioFormat RTAUDIO_SINT24
+    RtAudioFormat RTAUDIO_SINT32
+    RtAudioFormat RTAUDIO_FLOAT32
+    RtAudioFormat RTAUDIO_FLOAT64
+
+    ctypedef int (*RtAudioCallback)(void* outputBuffer,
+                                    void* inputBuffer,
+                                    unsigned int nFrames,
+                                    double streamTime,
+                                    RtAudioStreamStatus status,
+                                    void* userData)
+
+    ctypedef void (*RtAudioErrorCallback)(RtAudioError.Type _type,
+                                          const string& errortxt)
+
+    cdef cppclass cppRtAudio "RtAudio":
+        cppclass DeviceInfo:
+            bool probed
+            string name
+            unsigned int outputChannels
+            unsigned int inputChannels
+            unsigned int duplexChannels
+            bool isDefaultOutput
+            bool isDefaultInput
+            vector[unsigned int] sampleRates
+            unsigned int preferredSampleRate
+            RtAudioFormat nativeFormats
+
+        cppclass StreamOptions:
+            RtAudioStreamFlags flags
+            unsigned int numberOfBuffers
+            string streamName
+            int priority
+        cppclass StreamParameters:
+            unsigned int deviceId
+            unsigned int nChannels
+            unsigned int firstChannel
+
+        RtAudio() except +
+        # ~RtAudio() Destructors should not be listed
+        unsigned int getDeviceCount()
+        DeviceInfo getDeviceInfo(unsigned int device)
+        unsigned int getDefaultOutputDevice()
+        unsigned int getDefaultInputDevice()
+        void openStream(StreamParameters* outputParameters,
+                        StreamParameters* intputParameters,
+                        RtAudioFormat _format,
+                        unsigned int sampleRate,
+                        unsigned int* bufferFrames,
+                        RtAudioCallback callback,
+                        void* userData,
+                        void* StreamOptions,
+                        RtAudioErrorCallback) except +
+        void closeStream()
+        void startStream() except +
+        void stopStream() except +
+        void abortStream() except +
+        bool isStreamOpen()
+        bool isStreamRunning()
+        double getStreamTime()
+        void setStreamTime(double) except +
+        long getStreamLatency()
+        unsigned int getStreamSampleRate()
+        void showWarnings(bool value)
+    
+_formats_strkey = {
+        '8-bit integers': (RTAUDIO_SINT8,   1, np.int8),
+        '16-bit integers': (RTAUDIO_SINT16, 2, np.int16),
+        '24-bit integers': (RTAUDIO_SINT24, 3),
+        '32-bit integers': (RTAUDIO_SINT32, 4, np.int32),
+        '32-bit floats': (RTAUDIO_FLOAT32,  4, np.float32),
+        '64-bit floats': (RTAUDIO_FLOAT64,  8, np.float64),
+}
+_formats_rtkey = {
+        RTAUDIO_SINT8: ('8-bit integers', 1, cnp.NPY_INT8),
+        RTAUDIO_SINT16: ('16-bit integers',2, cnp.NPY_INT16),
+        RTAUDIO_SINT24: ('24-bit integers',3),
+        RTAUDIO_SINT32: ('32-bit integers',4, cnp.NPY_INT32),
+        RTAUDIO_FLOAT32: ('32-bit floats',  4, cnp.NPY_FLOAT32),
+        RTAUDIO_FLOAT64: ('64-bit floats',  8, cnp.NPY_FLOAT64),
+}
+
+def get_numpy_dtype_from_format_string(format_string):
+    return _formats_strkey[format_string][-1]
+def get_sampwidth_from_format_string(format_string):
+    return _formats_strkey[format_string][-2]
+
+cdef class _Stream:
+    cdef:
+        object pyCallback
+        unsigned int sampleSize
+        RtAudioFormat sampleformat
+        cppRtAudio.StreamParameters inputParams
+        cppRtAudio.StreamParameters outputParams
+        # These boolean values tell us whether the structs above here are
+        # initialized and contain valid data
+        bool hasInput
+        bool hasOutput
+        unsigned int bufferFrames
+
+
+# It took me quite a long time to fully understand Cython's idiosyncrasies
+# concerning C(++) callbacks, the GIL and passing Python objects as pointers
+# into C(++) functions. But finally, here it is!
+cdef object fromBufferToNPYNoCopy(
+                             cnp.NPY_TYPES buffer_format_type,
+                             void* buf,
+                             size_t nchannels,
+                             size_t nframes):
+    cdef cnp.npy_intp[2] dims = [nframes, nchannels]
+ 
+    # Interleaved data is C-style contiguous. Therefore, we can directly use
+    # SimpleNewFromData()
+    array = cnp.PyArray_SimpleNewFromData(2, &dims[0], buffer_format_type,
+            buf)
+
+    return array
+
+cdef void fromNPYToBuffer(cnp.ndarray arr,
+                          void* buf):
+    """
+    Copy a Python numpy array over to a buffer
+    No checks, just memcpy! Careful!
+    """
+    memcpy(buf, arr.data, arr.size*arr.itemsize)
+
+
+cdef int audioCallback(void* outputbuffer,
+                       void* inputbuffer,
+                       unsigned int nFrames,
+                       double streamTime,
+                       RtAudioStreamStatus status,
+                       void* userData) nogil:
+    """
+    Calls the Python callback function and converts data
+
+    """
+    cdef:
+        int rval = 0
+        cnp.NPY_TYPES npy_format
+
+    with gil:
+        if status == RTAUDIO_INPUT_OVERFLOW:
+            print('Input overflow.')
+            return 0
+        if status == RTAUDIO_OUTPUT_UNDERFLOW:
+            print('Output underflow.')
+            return 0
+        stream = <_Stream>(userData)
+        
+        # Obtain stream information
+        npy_input = None
+        npy_output = None
+        if stream.hasInput: 
+            try:
+                assert inputbuffer != NULL
+                npy_format = _formats_rtkey[stream.sampleformat][2]
+                npy_input = fromBufferToNPYNoCopy(
+                        npy_format,
+                        inputbuffer,
+                        stream.inputParams.nChannels,
+                        nFrames)
+
+            except Exception as e:
+                print('exception in cython callback: ', str(e))
+
+        if stream.hasOutput:
+            try:
+                assert outputbuffer != NULL
+                npy_format = _formats_rtkey[stream.sampleformat][2]
+                npy_output = fromBufferToNPYNoCopy(
+                        npy_format,
+                        outputbuffer,
+                        stream.outputParams.nChannels,
+                        nFrames)
+
+            except Exception as e:
+                print('exception in cython callback: ', str(e))
+        try:
+            rval = stream.pyCallback(npy_input,
+                                     npy_output,
+                                     nFrames,
+                                     streamTime)
+        except Exception as e:
+            print('Exception in Python callback: ', str(e))
+            return 1
+
+        return rval
+
+cdef void errorCallback(RtAudioError.Type _type,const string& errortxt) nogil:
+    pass
+
+
+
+
+cdef class RtAudio:
+    cdef:
+        cppRtAudio _rtaudio
+        _Stream _stream
+
+    def __cinit__(self):
+        self._stream = None
+
+    def __dealloc__(self):
+        if self._stream is not None:
+            print('Force closing stream')
+            self._rtaudio.closeStream()
+
+    cpdef unsigned int getDeviceCount(self):
+        return self._rtaudio.getDeviceCount()
+
+    cpdef unsigned int getDefaultOutputDevice(self):
+        return self._rtaudio.getDefaultOutputDevice()
+
+    cpdef unsigned int getDefaultInputDevice(self):
+        return self._rtaudio.getDefaultInputDevice()
+
+    def getDeviceInfo(self, unsigned int device):
+        """
+        Return device information of the current device
+        """
+        cdef cppRtAudio.DeviceInfo devinfo = self._rtaudio.getDeviceInfo(device)
+        sampleformats = []
+        nf = devinfo.nativeFormats
+        for format_ in [ RTAUDIO_SINT8, RTAUDIO_SINT16, RTAUDIO_SINT24,
+                RTAUDIO_SINT32, RTAUDIO_FLOAT32, RTAUDIO_FLOAT64]:
+            if nf & format_:
+                sampleformats.append(_formats_rtkey[format_][0])
+        return DeviceInfo(
+                index = device,
+                probed = devinfo.probed,
+                name = devinfo.name.decode('utf-8'),
+                outputchannels = devinfo.outputChannels,
+                inputchannels = devinfo.inputChannels,
+                duplexchannels = devinfo.duplexChannels,
+                samplerates = devinfo.sampleRates,
+                sampleformats = sampleformats,
+                prefsamplerate = devinfo.preferredSampleRate)
+
+    @cython.nonecheck(True)
+    def openStream(self,object outputParams,
+                        object inputParams,
+                        str sampleformat,
+                        unsigned int sampleRate,
+                        unsigned int bufferFrames,
+                        object pyCallback,
+                        object options = None,
+                        object pyErrorCallback = None):
+        """
+        Opening a stream with specified parameters
+
+        Args:
+            outputParams: dictionary of stream outputParameters, set to None
+            if no outputPararms are specified
+            inputParams: dictionary of stream inputParameters, set to None
+            if no inputPararms are specified
+            sampleRate: desired sample rate.
+            bufferFrames: the amount of frames in a callback buffer
+            callback: callback to call. Note: this callback is called on a
+            different thread!
+            options: A dictionary of optional additional stream options
+            errorCallback: client-defined function that will be invoked when an
+            error has occured.
+
+        Returns: None
+        """
+        if self._stream is not None:
+            raise RuntimeError('Stream is already opened.')
+
+        cdef cppRtAudio.StreamParameters *rtOutputParams_ptr = NULL
+        cdef cppRtAudio.StreamParameters *rtInputParams_ptr = NULL
+
+        cdef cppRtAudio.StreamOptions streamoptions
+        streamoptions.flags = RTAUDIO_HOG_DEVICE
+        streamoptions.numberOfBuffers = 4
+
+        self._stream = _Stream()
+        self._stream.pyCallback = pyCallback
+        self._stream.sampleformat = _formats_strkey[sampleformat][0]
+
+        if outputParams is not None:
+            rtOutputParams_ptr = &self._stream.outputParams
+            rtOutputParams_ptr.deviceId = outputParams['deviceid']
+            rtOutputParams_ptr.nChannels = outputParams['nchannels']
+            rtOutputParams_ptr.firstChannel = outputParams['firstchannel']
+            self._stream.hasOutput = True
+
+        if inputParams is not None:
+            rtInputParams_ptr = &self._stream.inputParams
+            rtInputParams_ptr.deviceId = inputParams['deviceid']
+            rtInputParams_ptr.nChannels = inputParams['nchannels']
+            rtInputParams_ptr.firstChannel = inputParams['firstchannel']
+            self._stream.hasInput = True
+
+        try:
+            self._stream.bufferFrames = bufferFrames
+            self._rtaudio.openStream(rtOutputParams_ptr,
+                                     rtInputParams_ptr,
+                                     _formats_strkey[sampleformat][0],
+                                     sampleRate,
+                                     &self._stream.bufferFrames,
+                                     audioCallback,
+                                     <void*> self._stream,
+                                     &streamoptions, # Stream options
+                                     errorCallback  # Error callback
+                                     ) 
+        except Exception as e:
+            print('Exception occured in stream opening: ', str(e))
+            self._stream = None
+            raise
+
+        return self._stream.bufferFrames
+
+    def startStream(self):
+        self._rtaudio.startStream()
+
+    def stopStream(self):
+        if not self._stream:
+            raise RuntimeError('Stream is not opened')
+        self._rtaudio.stopStream()
+
+    def closeStream(self):
+        if not self._stream:
+            raise RuntimeError('Stream is not opened')
+        # Closing stream
+        self._rtaudio.closeStream()
+        self._stream = None
+
+    def abortStream(self):
+        if not self._stream:
+            raise RuntimeError('Stream is not opened')
+        self._rtaudio.abortStream()
+
+    def isStreamOpen(self):
+        return self._rtaudio.isStreamOpen()
+
+    def isStreamRunning(self):
+        return self._rtaudio.isStreamRunning()
+
+    def getStreamTime(self):
+        return self._rtaudio.getStreamTime()
+    
+    def setStreamTime(self, double time):
+        return self._rtaudio.setStreamTime(time)
+
+    

lasp/filter/__init__.py

@@ -1,4 +1,4 @@
+from .soundpressureweighting import *
+from .filterbank_design import *
+from .fir_design import *
 
-from .freqweighting_fir import A, C
-
-__all__ = ['A', 'Z']

lasp/filter/bandpass_fir.py

@@ -1,234 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""!
-Author: J.A. de Jong - ASCEE
-
-Description: FIR filter design for octave bands from 16Hz to 16 kHz for a
-sampling frequency of 48 kHz, FIR filter design for one-third octave bands
-
-See test/octave_fir_test.py for a testing
-
-"""
-from .fir_design import bandpass_fir_design, freqResponse as frsp
-import numpy as np
-__all__ = ['OctaveBankDesigner', 'ThirdOctaveBankDesigner']
-
-
-class FilterBankDesigner:
-    """
-    A class responsible for designing FIR filters
-    """
-    G = 10**(3/10)
-    fr = 1000.
-    L = 256  # Filter order
-    fs = 48000.  # Sampling frequency
-
-    def fm(self, x):
-        """
-        Returns the exact midband frequency of the bandpass filter
-
-        Args:
-            x: Midband designator
-        """
-        # Exact midband frequency
-        return self.G**(x/self.b)*self.fr
-
-    def fl(self, x):
-        """
-        Returns the exact cut-on frequency of the bandpass filter
-
-        Args:
-            x: Midband designator
-        """
-        return self.fm(x)*self.G**(-1/(2*self.b))
-
-    def fu(self, x):
-        """
-        Returns the exact cut-off frequency of the bandpass filter
-
-        Args:
-            x: Midband designator
-        """
-        return self.fm(x)*self.G**(1/(2*self.b))
-
-    def createFilter(self, fs, x):
-        """
-        Create a FIR filter for band designator b and sampling frequency fs.
-        Decimation should be obtained from decimation() method.
-
-        Returns:
-            filter: 1D ndarray with FIR filter coefficients
-        """
-        assert np.isclose(fs, self.fs), "Invalid sampling frequency"
-        fd = fs / np.prod(self.decimation(x))
-
-        # For designing the filter, the lower and upper frequencies need to be
-        # slightly adjusted to fall within the limits for a class 1 filter.
-        fl = self.fl(x)*self.fac_l(x)
-        fu = self.fu(x)*self.fac_u(x)
-
-        return bandpass_fir_design(self.L, fd, fl, fu)
-
-    def freqResponse(self, fs, x, freq):
-        """
-        Compute the frequency response for a certain filter
-
-        Args:
-            fs: Sampling frequency [Hz]
-            x: Midband designator
-        """
-        fir = self.createFilter(fs, x)
-
-        # Decimated sampling frequency [Hz]
-        fd = fs / np.prod(self.decimation(x))
-
-        return frsp(fd, freq, fir)
-
-    def nominal_txt_tox(self, nom_txt):
-        """
-        Returns the x-value corresponding to a certain nominal txt: '1k' -> 0
-        """
-        for x in self.xs:
-            if self.nominal_txt(x) == nom_txt:
-                return x
-        raise ValueError(
-            f'Could not find an x-value corresponding to {nom_txt}.')
-
-
-class OctaveBankDesigner(FilterBankDesigner):
-    """
-    Octave band filter designer
-    """
-
-    def __init__(self):
-        pass
-
-    @property
-    def b(self):
-        # Band division, 1 for octave bands
-        return 1
-
-    @property
-    def xs(self):
-        return list(range(-6, 5))
-
-    def nominal_txt(self, x):
-        # Text corresponding to the nominal frequency
-        nominals = {4: '16k',
-                    3: '8k',
-                    2: '4k',
-                    1: '2k',
-                    0: '1k',
-                    -1: '500',
-                    -2: '250',
-                    -3: '125',
-                    -4: '63',
-                    -5: '31.5',
-                    -6: '16'}
-        return nominals[x]
-
-    def fac_l(self, x):
-        """
-        Factor with which to multiply the cut-on frequency of the FIR filter
-        """
-        if x == 4:
-            return .995
-        elif x in (3, 1):
-            return .99
-        elif x in(-6, -4, -2, 2, 0):
-            return .98
-        else:
-            return .96
-
-    def fac_u(self, x):
-        """
-        Factor with which to multiply the cut-off frequency of the FIR filter
-        """
-        if x == 4:
-            return 1.004
-        elif x in (3, 1):
-            return 1.006
-        elif x in (-6, -4, -2, 2, 0):
-            return 1.01
-        else:
-            return 1.02
-
-    def decimation(self, x):
-        """
-        Required decimation for each filter
-        """
-        if x > 1:
-            return [1]
-        elif x > -2:
-            return [4]
-        elif x > -4:
-            return [4, 4]
-        elif x > -6:
-            return [4, 4, 4]
-        elif x == -6:
-            return [4, 4, 4, 4]
-        assert False, 'Overlooked decimation'
-
-
-class ThirdOctaveBankDesigner(FilterBankDesigner):
-
-    def __init__(self):
-        self.xs = list(range(-16, 14))
-        # Text corresponding to the nominal frequency
-        self._nominal_txt = ['25', '31.5', '40',
-                            '50', '63', '80',
-                            '100', '125', '160',
-                            '200', '250', '315',
-                            '400', '500', '630',
-                            '800', '1k', '1.25k',
-                            '1.6k', '2k', '2.5k',
-                            '3.15k', '4k', '5k',
-                            '6.3k', '8k', '10k',
-                            '12.5k', '16k', '20k']
-
-        assert len(self.xs) == len(self._nominal_txt)
-
-    @property
-    def b(self):
-        # Band division factor, 3 for one-third octave bands
-        return 3
-
-    def nominal_txt(self, x):
-        # Put the nominal frequencies in a dictionary for easy access with
-        # x as the key.
-        index = x - self.xs[0]
-        return self._nominal_txt[index]
-
-    @staticmethod
-    def decimation(x):
-        if x > 5:
-            return [1]
-        elif x > -1:
-            return [4]
-        elif x > -7:
-            return [4, 4]
-        elif x > -13:
-            return [4, 4, 4]
-        elif x > -17:
-            return [4, 4, 4, 4]
-        assert False, 'Bug: overlooked decimation'
-
-    @staticmethod
-    def fac_l(x):
-        if x in (-13, -7, -1, 5, 11, 12, 13):
-            return .995
-        elif x in (-12, -6, 0, 6):
-            return .98
-        else:
-            return .99
-
-    @staticmethod
-    def fac_u(x):
-        if x in (-14, -13, -8, -7, -1, -2, 3, 4, 5, 10, 11, 12):
-            return 1.005
-        elif x in (12, 13):
-            return 1.003
-        elif x in (12, -6, 0):
-            return 1.015
-        else:
-            return 1.01

lasp/filter/bandpass_limits.py

@@ -1,99 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""!
-Author: J.A. de Jong - ASCEE
-
-Description: Limit lines for class 1 octave band filter limits according to
-the ICS 17.140.50 standard.
-"""
-__all__ = ['G', 'fr', 'third_octave_band_limits', 'octave_band_limits']
-import numpy as np
-
-# Reference frequency
-fr = 1000.
-G = 10**(3/10)
-
-
-def third_octave_band_limits(x):
-    """
-    Returns the class 1 third octave band filter limits for filter designator
-    x.
-
-    Args:
-        x: Filter offset power from the reference frequency of 1000 Hz.
-
-    Returns:
-        freq, ulim, llim: Tuple of Numpy arrays containing the frequencyies,
-        upper and lower limits of the filter.
-    """
-    b = 3
-
-    fm = G**(x/b)*fr
-    plusinf = 20
-    f_ratio_pos = [1., 1.02667, 1.05575, 1.08746, 1.12202, 1.12202,
-                   1.29437, 1.88173, 3.05365, 5.39195, plusinf]
-
-    f_ratio_neg = [0.97402, 0.94719, 0.91958, 0.89125, 0.89125,
-                   0.77257, 0.53143, 0.32748, 0.18546, 1/plusinf]
-    f_ratio_neg.reverse()
-
-    f_ratio = f_ratio_neg + f_ratio_pos
-
-    mininf = -1e300
-
-    upper_limits_pos = [.3]*5 + [-2, -17.5, -42, -61, -70, -70]
-    upper_limits_neg = upper_limits_pos[:]
-    upper_limits_neg.reverse()
-    upper_limits = np.array(upper_limits_neg[:-1] + upper_limits_pos)
-
-    lower_limits_pos = [-.3, -.4, -.6, -1.3, -5, -5, mininf, mininf,
-                        mininf, mininf, mininf]
-    lower_limits_neg = lower_limits_pos[:]
-    lower_limits_neg.reverse()
-    lower_limits = np.array(lower_limits_neg[:-1] + lower_limits_pos)
-
-    freqs = fm*np.array(f_ratio)
-
-    return freqs, upper_limits, lower_limits
-
-
-def octave_band_limits(x):
-
-    b = 1
-
-    # Exact midband frequency
-    fm = G**(x/b)*fr
-
-    G_power_values_pos = [0, 1/8, 1/4, 3/8, 1/2, 1/2, 1, 2, 3, 4]
-    G_power_values_neg = [-i for i in G_power_values_pos]
-    G_power_values_neg.reverse()
-    G_power_values = G_power_values_neg[:-1] + G_power_values_pos
-
-    mininf = -1e300
-
-    lower_limits_pos = [-0.3, -0.4, -0.6, -1.3, -5.0, -5.0] + 4*[mininf]
-    lower_limits_neg = lower_limits_pos[:]
-    lower_limits_neg.reverse()
-    lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos)
-
-    upper_limits_pos = [0.3]*5 + [-2, -17.5, -42, -61, -70]
-    upper_limits_neg = upper_limits_pos[:]
-    upper_limits_neg.reverse()
-    upper_limits = np.asarray(upper_limits_neg[:-1] + upper_limits_pos)
-
-    freqs = fm*G**np.asarray(G_power_values)
-
-    return freqs, upper_limits, lower_limits
-
-
-if __name__ == '__main__':
-
-    from asceefigs.plot import close, Figure
-    close('all')
-    freqs, upper_limits, lower_limits = octave_band_limits(0)
-
-    f = Figure()
-    f.semilogx(freqs, lower_limits)
-    f.semilogx(freqs, upper_limits)
-
-    f.ylim(-80, 1)

lasp/filter/filterbank_design.py

@@ -0,0 +1,445 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE
+
+Description: FIR filter design for octave bands from 16Hz to 16 kHz for a
+sampling frequency of 48 kHz, filter design for one-third octave bands.
+Resulting filters are supposed to be standard compliant.
+
+See test/octave_fir_test.py for a testing
+
+"""
+from .fir_design import bandpass_fir_design, freqResponse as firFreqResponse
+import numpy as np
+
+# For designing second-order sections
+from scipy.signal import butter
+
+
+__all__ = ['OctaveBankDesigner', 'ThirdOctaveBankDesigner']
+
+
+class FilterBankDesigner:
+    """A class responsible for designing FIR filters."""
+
+    def __init__(self, fs):
+        """Initialize a filter bank designer.
+
+        Args:
+            fs: Sampling frequency [Hz]
+        """
+        # Default FIR filter length
+        firFilterLength = 256  # Filter order
+        self.fs = fs
+
+        # Constant G, according to standard
+        self.G = 10**(3/10)
+
+        # Reference frequency for all filter banks
+        self.fr = 1000.
+
+    def testStandardCompliance(self, x, freq, h_dB, filter_class=0):
+        """Test whether the filter with given frequency response is compliant 
+        with the standard.
+
+        Args:
+            x: Band designator
+            freq: Array of frequencies to test for. Note: *Should be fine
+            enough to follow response!*
+            h_dB: Filter frequency response in *deciBell*
+            filter_class: Filter class to test for
+
+        Returns:
+            True if filter is norm-compliant, False if not
+        
+        """
+        # Skip zero-frequency
+        if np.isclose(freq[0], 0):
+            freq = freq[1:]
+            h_dB = h_dB[1:]
+        freqlim, llim, ulim = self.band_limits(x, filter_class)
+
+        # Interpolate limites to frequency array as given
+        llim_full = np.interp(freq, freqlim, llim, left=-np.inf, right=-np.inf)
+        ulim_full = np.interp(freq, freqlim, ulim, left=ulim[0], right=ulim[-1])
+
+        return bool(np.all(llim_full <= h_dB) and
+                    np.all(ulim_full >= h_dB))
+
+    def fm(self, x):
+        """Returns the exact midband frequency of the bandpass filter.
+
+        Args:
+            x: Midband designator
+        """
+        # Exact midband frequency
+        return self.G**(x/self.b)*self.fr
+
+    def fl(self, x):
+        """Returns the exact cut-on frequency of the bandpass filter.
+
+        Args:
+            x: Midband designator
+        """
+        return self.fm(x)*self.G**(-1/(2*self.b))
+
+    def fu(self, x):
+        """Returns the exact cut-off frequency of the bandpass filter.
+
+        Args:
+            x: Midband designator
+        """
+        return self.fm(x)*self.G**(1/(2*self.b))
+
+    def createFirFilter(self, x):
+        """Create a FIR filter for band designator b and sampling frequency fs.
+        firdecimation should be obtained from firdecimation() method.
+
+        Returns:
+            filter: 1D ndarray with FIR filter coefficients
+        """
+        assert np.isclose(fs, self.fs), "Invalid sampling frequency"
+        fd = fs / np.prod(self.firDecimation(x))
+
+        # For designing the filter, the lower and upper frequencies need to be
+        # slightly adjusted to fall within the limits for a class 1 filter.
+        fl = self.fl(x)*self.firFac_l(x)
+        fu = self.fu(x)*self.firFac_u(x)
+
+        return bandpass_fir_design(self.firFilterLength, fd, fl, fu)
+
+    def createSOSFilter(self, x: int):
+        """Create a Second Order Section filter (cascaded BiQuad's) for the
+        given sample rate and band designator.
+
+        Args:
+            x: Band designator
+        """
+
+        SOS_ORDER = 5
+
+        fs = self.fs
+        fl = self.fl(x)*self.sosFac_l(x)
+        fu = self.fu(x)*self.sosFac_u(x)
+
+        fnyq = fs/2
+
+        # Normalized upper and lower frequencies of the bandpass
+        fl_n = fl/fnyq
+        fu_n = fu/fnyq
+
+        return butter(SOS_ORDER, [fl_n, fu_n], output='sos', btype='band')
+
+    def firFreqResponse(self, x, freq):
+        """Compute the frequency response for a certain filter.
+
+        Args:
+            x: Midband designator
+            freq: Array of frequencies to evaluate on
+
+        Returns:
+            h: Linear filter transfer function [-]
+        """
+        fir = self.createFirFilter(fs, x)
+
+        # Decimated sampling frequency [Hz]
+        fd = fs / np.prod(self.firdecimation(x))
+
+        return firFreqResponse(fd, freq, fir)
+
+    def nominal_txt_tox(self, nom_txt: str):
+        """Returns the x-value corresponding to a certain nominal txt: '1k' ->
+        0.
+
+        Args:
+            nom_txt: Text-representation of midband frequency
+        """
+        for x in self.xs:
+            if self.nominal_txt(x) == nom_txt:
+                return x
+        raise ValueError(
+            f'Could not find an x-value corresponding to {nom_txt}.')
+
+
+class OctaveBankDesigner(FilterBankDesigner):
+    """Octave band filter designer."""
+
+    def __init__(self, fs):
+        super().__init__(fs)
+
+    @property
+    def b(self):
+        # Band division, 1 for octave bands
+        return 1
+
+    @property
+    def xs(self):
+        """All possible band designators for an octave band filter."""
+        return list(range(-6, 5))
+
+    def band_limits(self, x, filter_class=0):
+        """Returns the octave band filter limits for filter designator x.
+
+        Args:
+            x: Filter offset power from the reference frequency of 1000 Hz.
+            filter_class: Either 0 or 1, defines the tolerances on the frequency
+            response
+
+        Returns:
+            freq, llim, ulim: Tuple of Numpy arrays containing the frequencies of
+            the corner points of the filter frequency response limits, lower limits
+            in *deciBell*, upper limits in *deciBell*, respectively.
+        """
+        b = 1
+
+        # Exact midband frequency
+        fm = self.G**(x/self.b)*self.fr
+
+        G_power_values_pos = [0, 1/8, 1/4, 3/8, 1/2, 1/2, 1, 2, 3, 4]
+        G_power_values_neg = [-i for i in G_power_values_pos]
+        G_power_values_neg.reverse()
+        G_power_values = G_power_values_neg[:-1] + G_power_values_pos
+
+        mininf = -1e300
+
+        if filter_class == 1:
+            lower_limits_pos = [-0.3, -0.4, -0.6, -1.3, -5.0, -5.0] + 4*[mininf]
+        elif filter_class == 0:
+            lower_limits_pos = [-0.15, -0.2, -0.4, -1.1, -4.5, -4.5] + 4*[mininf]
+        lower_limits_neg = lower_limits_pos[:]
+        lower_limits_neg.reverse()
+        lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos)
+
+        if filter_class == 1:
+            upper_limits_pos = [0.3]*5 + [-2, -17.5, -42, -61, -70]
+        if filter_class == 0:
+            upper_limits_pos = [0.15]*5 + [-2.3, -18, -42.5, -62, -75]
+        upper_limits_neg = upper_limits_pos[:]
+        upper_limits_neg.reverse()
+        upper_limits = np.asarray(upper_limits_neg[:-1] + upper_limits_pos)
+
+        freqs = fm*self.G**np.asarray(G_power_values)
+
+        return freqs, lower_limits, upper_limits
+
+
+    def nominal_txt(self, x):
+        """Returns textual repressentation of corresponding to the nominal
+        frequency."""
+        nominals = {4: '16k',
+                    3: '8k',
+                    2: '4k',
+                    1: '2k',
+                    0: '1k',
+                    -1: '500',
+                    -2: '250',
+                    -3: '125',
+                    -4: '63',
+                    -5: '31.5',
+                    -6: '16'}
+        assert len(nominals) == len(self.xs)
+        return nominals[x]
+
+    def firFac_l(self, x):
+        """Factor with which to multiply the cut-on frequency of the FIR
+        filter."""
+        assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
+        if x == 4:
+            return .995
+        elif x in (3, 1):
+            return .99
+        elif x in(-6, -4, -2, 2, 0):
+            return .98
+        else:
+            return .96
+
+    def firFac_u(self, x):
+        """Factor with which to multiply the cut-off frequency of the FIR
+        filter."""
+        assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
+        if x == 4:
+            return 1.004
+        elif x in (3, 1):
+            return 1.006
+        elif x in (-6, -4, -2, 2, 0):
+            return 1.01
+        else:
+            return 1.02
+
+    def firDecimation(self, x):
+        """Required firdecimation for each filter."""
+        assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
+        if x > 1:
+            return [1]
+        elif x > -2:
+            return [4]
+        elif x > -4:
+            return [4, 4]
+        elif x > -6:
+            return [4, 4, 4]
+        elif x == -6:
+            return [4, 4, 4, 4]
+        assert False, 'Overlooked firdecimation'
+
+    def sosFac_l(self, x):
+        """Left side percentage of change in cut-on frequency for designing the
+        filter, for OCTAVE band filter.
+
+        Args:
+            x: Filter band designator
+        """
+        # Idea: correct for frequency warping:
+        if int(self.fs) in [48000, 96000]:
+            return 1.0
+        else:
+            raise ValueError('Unimplemented sampling frequency for SOS'
+                             'filter design')
+
+    def sosFac_u(self, x):
+        """Right side percentage of change in cut-on frequency for designing
+        the filter.
+
+        Args:
+            x: Filter band designator
+        """
+        if int(self.fs) in [48000, 96000]:
+            return 1.0
+        else:
+            raise ValueError('Unimplemented sampling frequency for SOS'
+                             'filter design')
+
+
+class ThirdOctaveBankDesigner(FilterBankDesigner):
+
+    def __init__(self, fs):
+        super().__init__(fs)
+        self.xs = list(range(-16, 14))
+        # Text corresponding to the nominal frequency
+        self._nominal_txt = ['25', '31.5', '40',
+                             '50', '63', '80',
+                             '100', '125', '160',
+                             '200', '250', '315',
+                             '400', '500', '630',
+                             '800', '1k', '1.25k',
+                             '1.6k', '2k', '2.5k',
+                             '3.15k', '4k', '5k',
+                             '6.3k', '8k', '10k',
+                             '12.5k', '16k', '20k']
+
+        assert len(self.xs) == len(self._nominal_txt)
+
+    @property
+    def b(self):
+        # Band division factor, 3 for one-third octave bands
+        return 3
+
+    def nominal_txt(self, x):
+        # Put the nominal frequencies in a dictionary for easy access with
+        # x as the key.
+        index = x - self.xs[0]
+        return self._nominal_txt[index]
+
+    def band_limits(self, x, filter_class=0):
+        """Returns the third octave band filter limits for filter designator x.
+
+        Args:
+            x: Filter offset power from the reference frequency of 1000 Hz.
+            filter_class: Either 0 or 1, defines the tolerances on the frequency
+            response
+
+        Returns:
+            freq, llim, ulim: Tuple of Numpy arrays containing the frequencies of
+            the corner points of the filter frequency response limits, lower limits
+            in *deciBell*, upper limits in *deciBell*, respectively.
+        """
+
+        fm = self.G**(x/self.b)*self.fr
+        plusinf = 20
+        f_ratio_pos = [1., 1.02667, 1.05575, 1.08746, 1.12202, 1.12202,
+                       1.29437, 1.88173, 3.05365, 5.39195, plusinf]
+
+        f_ratio_neg = [0.97402, 0.94719, 0.91958, 0.89125, 0.89125,
+                       0.77257, 0.53143, 0.32748, 0.18546, 1/plusinf]
+        f_ratio_neg.reverse()
+
+        f_ratio = f_ratio_neg + f_ratio_pos
+
+        mininf = -1e300
+
+        if filter_class == 1:
+            upper_limits_pos = [.3]*5 + [-2, -17.5, -42, -61, -70, -70]
+        elif filter_class == 0:
+            upper_limits_pos = [.15]*5 + [-2.3, -18, -42.5, -62, -75, -75]
+        else:
+            raise ValueError('Filter class should either be 0 or 1')
+
+        upper_limits_neg = upper_limits_pos[:]
+        upper_limits_neg.reverse()
+        upper_limits = np.array(upper_limits_neg[:-1] + upper_limits_pos)
+
+        if filter_class == 1:
+            lower_limits_pos = [-.3, -.4, -.6, -1.3, -5, -5, mininf, mininf,
+                                mininf, mininf, mininf]
+        elif filter_class == 0:
+            lower_limits_pos = [-.15, -.2, -.4, -1.1, -4.5, -4.5, mininf, mininf,
+                                mininf, mininf, mininf]
+
+        lower_limits_neg = lower_limits_pos[:]
+        lower_limits_neg.reverse()
+        lower_limits = np.array(lower_limits_neg[:-1] + lower_limits_pos)
+
+        freqs = fm*np.array(f_ratio)
+
+        return freqs, lower_limits, upper_limits
+
+    def firDecimation(self, x):
+        assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
+        if x > 5:
+            return [1]
+        elif x > -1:
+            return [4]
+        elif x > -7:
+            return [4, 4]
+        elif x > -13:
+            return [4, 4, 4]
+        elif x > -17:
+            return [4, 4, 4, 4]
+        assert False, 'Bug: overlooked firdecimation'
+
+    def firFac_l(self, x):
+        if x in (-13, -7, -1, 5, 11, 12, 13):
+            return .995
+        elif x in (-12, -6, 0, 6):
+            return .98
+        else:
+            return .99
+
+    def firFac_u(self, x):
+        if x in (-14, -13, -8, -7, -1, -2, 3, 4, 5, 10, 11, 12):
+            return 1.005
+        elif x in (12, 13):
+            return 1.003
+        elif x in (12, -6, 0):
+            return 1.015
+        else:
+            return 1.01
+
+    def sosFac_l(self, x):
+        """Left side percentage of change in cut-on frequency for designing the
+        filter."""
+        # Idea: correct for frequency warping:
+        if np.isclose(self.fs, 48000):
+            return 1.00
+        else:
+            raise ValueError('Unimplemented sampling frequency for SOS'
+                             'filter design')
+
+    def sosFac_u(self, x):
+        """Right side percentage of change in cut-on frequency for designing
+        the filter."""
+        if np.isclose(self.fs, 48000):
+            return 1
+        else:
+            raise ValueError('Unimplemented sampling frequency for SOS'
+                             'filter design')

lasp/filter/soundpressureweighting.py

@@ -3,7 +3,8 @@
 """!
 Author: J.A. de Jong - ASCEE
 
-Description: Filter design for frequency weightings A and C.
+Description: Filter design for frequency weighting curves (i.e. A and C 
+weighting)
 """
 from .fir_design import freqResponse, arbitrary_fir_design
 import numpy as np
@@ -105,8 +106,7 @@ def C_fir_design():
 
 
 def show_Afir():
-    from asceefigs.plot import close, Figure
-    close('all')
+    from asceefig.plot import Figure
 
     fs = 48000.
     freq_design = np.linspace(0, 17e3, 3000)
@@ -132,8 +132,7 @@ def show_Afir():
 
 
 def show_Cfir():
-    from asceefigs.plot import close, Figure
-    close('all')
+    from asceefig.plot import Figure
 
     fs = 48000.
     freq_design = np.linspace(0, 17e3, 3000)

lasp/lasp_avstream.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python3.6
 # -*- coding: utf-8 -*-
 """
 Description: Read data from image stream and record sound at the same time
@@ -6,36 +6,104 @@ Description: Read data from image stream and record sound at the same time
 import cv2 as cv
 from .lasp_atomic import Atomic
 from threading import Thread, Condition, Lock
+
 import time
-from .device import DAQDevice, roga_plugndaq
-import numpy as np
+from .device import (RtAudio, DeviceInfo, DAQConfiguration,
+                     get_numpy_dtype_from_format_string,
+                     get_sampwidth_from_format_string)
+
 __all__ = ['AvType', 'AvStream']
 
 video_x, video_y = 640, 480
-dtype, sampwidth = 'int16', 2
 
 
 class AvType:
-    video = 0
-    audio = 1
+    """
+    Specificying the type of data, for adding and removing callbacks from the
+    stream.
+    """
+    audio_input = 0
+    audio_output = 1
+    video = 2
 
 
 class AvStream:
-    def __init__(self, daqconfig=roga_plugndaq, video=None):
+    """
+    Audio and video data stream, to which callbacks can be added for processing
+    the data.
+    """
+
+    def __init__(self,
+                 device: DeviceInfo,
+                 avtype: AvType,
+                 daqconfig: DAQConfiguration,
+                 video=None):
+        """
+        Open a stream for audio in/output and video input. For audio output, by
+        default all available channels are opened for outputting data.
+
+        Args:
+            device: DeviceInfo for the audio device
+            avtype: Type of stream. Input, output or duplex
+
+            daqconfig: DAQConfiguration instance. If duplex mode flag is set,
+            please make sure that no output_device is None, as in that case the
+            output config will be taken from the input device.
+            video:
+        """
 
         self.daqconfig = daqconfig
+        self._device = device
+        self.avtype = avtype
+        self.duplex_mode = daqconfig.duplex_mode
+
+        # Determine highest input channel number
+        channelconfigs = daqconfig.en_input_channels
+
+        self.sensitivity = self.daqconfig.getSensitivities()
+
+        rtaudio_inputparams = None
+        rtaudio_outputparams = None
+
+        self.nframes_per_block = 2048
+
+        if self.duplex_mode or avtype == AvType.audio_output:
+            rtaudio_outputparams = {'deviceid': device.index,
+            # TODO: Add option to specify the number of output channels to use
+                                    'nchannels': 1,  #device.outputchannels,
+                                    'firstchannel': 0}
+            self.sampleformat = daqconfig.en_output_sample_format
+            self.samplerate = int(daqconfig.en_output_rate)
+
+        if avtype == AvType.audio_input or self.duplex_mode:
+            for i, channelconfig in enumerate(channelconfigs):
+                if channelconfig.channel_enabled:
+                    self.nchannels = i+1
+            rtaudio_inputparams = {'deviceid': device.index,
+                                   'nchannels': self.nchannels,
+                                   'firstchannel': 0}
+
+            # Here, we override the sample format in case of duplex mode.
+            self.sampleformat = daqconfig.en_input_sample_format
+            self.samplerate = int(daqconfig.en_input_rate)
+
         try:
-            daq = DAQDevice(daqconfig)
-            self.nchannels = len(daq.channels_en)
-            self.samplerate = daq.input_rate
-            self.blocksize = daq.blocksize
-            self.sensitivity = np.asarray(daqconfig.input_sensitivity)[
-                                          daq.channels_en]
+            self._rtaudio = RtAudio()
+            self.blocksize = self._rtaudio.openStream(
+                    rtaudio_outputparams,  # Outputparams
+                    rtaudio_inputparams,   # Inputparams
+                    self.sampleformat,     # Sampleformat
+                    self.samplerate,
+                    self.nframes_per_block, # Buffer size in frames
+                    self._audioCallback)
+
         except Exception as e:
             raise RuntimeError(f'Could not initialize DAQ device: {str(e)}')
 
-        self.video_x, self.video_y = video_x, video_y
-        self.dtype, self.sampwidth = dtype, sampwidth
+        self.numpy_dtype = get_numpy_dtype_from_format_string(
+            self.sampleformat)
+        self.sampwidth = get_sampwidth_from_format_string(
+            self.sampleformat)
 
         self._aframectr = Atomic(0)
         self._vframectr = Atomic(0)
@@ -47,28 +115,40 @@ class AvStream:
 
         self._video = video
         self._video_started = Atomic(False)
-        self._callbacks = []
-        self._audiothread = None
+        self._callbacks = {
+            AvType.audio_input: [],
+            AvType.audio_output: [],
+            AvType.video: []
+        }
         self._videothread = None
 
+    def close(self):
+        self._rtaudio.closeStream()
+
     def nCallbacks(self):
         """
         Returns the current number of installed callbacks
         """
-        return len(self._callbacks)
+        return len(self._callbacks[AvType.audio_input]) + \
+               len(self._callbacks[AvType.audio_output]) + \
+               len(self._callbacks[AvType.video]) 
 
-    def addCallback(self, cb):
+    def addCallback(self, cb: callable, cbtype: AvType):
         """
         Add as stream callback to the list of callbacks
         """
         with self._callbacklock:
-            if cb not in self._callbacks:
-                self._callbacks.append(cb)
+            outputcallbacks = self._callbacks[AvType.audio_output]
+            if cbtype == AvType.audio_output and len(outputcallbacks) > 0:
+                raise RuntimeError('Only one audio output callback can be allowed')
 
-    def removeCallback(self, cb):
+            if cb not in self._callbacks[cbtype]:
+                self._callbacks[cbtype].append(cb)
+
+    def removeCallback(self, cb, cbtype: AvType):
         with self._callbacklock:
-            if cb in self._callbacks:
-                self._callbacks.remove(cb)
+            if cb in self._callbacks[cbtype]:
+                self._callbacks[cbtype].remove(cb)
 
     def start(self):
         """
@@ -79,31 +159,15 @@ class AvStream:
         if self._running:
             raise RuntimeError('Stream already started')
 
-        assert self._audiothread is None
         assert self._videothread is None
 
         self._running <<= True
-        self._audiothread = Thread(target=self._audioThread)
         if self._video is not None:
             self._videothread = Thread(target=self._videoThread)
             self._videothread.start()
         else:
             self._video_started <<= True
-        self._audiothread.start()
-
-    def _audioThread(self):
-        # Raw stream to allow for in24 packed data type
-        try:
-            daq = DAQDevice(self.daqconfig)
-            # Get a single block first and do not process it. This one often
-            # contains quite some rubbish.
-            data = daq.read()
-            while self._running:
-                # print('read data...')
-                data = daq.read()
-                self._audioCallback(data)
-        except RuntimeError as e:
-            print(f'Runtime error occured during audio capture: {str(e)}')
+        self._rtaudio.startStream()
 
     def _videoThread(self):
         cap = cv.VideoCapture(self._video)
@@ -118,42 +182,50 @@ class AvStream:
                 if vframectr == 0:
                     self._video_started <<= True
                 with self._callbacklock:
-                    for cb in self._callbacks:
-                        cb(AvType.video, frame, self._aframectr(), vframectr)
+                    for cb in self._callbacks[AvType.video]:
+                        cb(frame, vframectr)
                 vframectr += 1
                 self._vframectr += 1
             else:
-
+                loopctr += 1
                 if loopctr == 10:
                     print('Error: no video capture!')
             time.sleep(0.2)
-            loopctr += 1
 
         cap.release()
         print('stopped videothread')
 
-    def _audioCallback(self, indata):
-        """This is called (from a separate thread) for each audio block."""
-        if not self._video_started:
-            return
-
-        with self._callbacklock:
-            for cb in self._callbacks:
-                cb(AvType.audio, indata, self._aframectr(), self._vframectr())
+    def _audioCallback(self, indata, outdata, nframes, streamtime):
+        """
+        This is called (from a separate thread) for each audio block.
+        """
         self._aframectr += 1
+        with self._callbacklock:
+            for cb in self._callbacks[AvType.audio_input]:
+                try:
+                    cb(indata, outdata, self._aframectr())
+                except Exception as e:
+                    print(e)
+                    return 1
+            for cb in self._callbacks[AvType.audio_output]:
+                try:
+                    cb(indata, outdata, self._aframectr())
+                except Exception as e:
+                    print(e)
+                    return 1
+        return 0 if self._running else 1
 
     def stop(self):
         self._running <<= False
         with self._running_cond:
             self._running_cond.notify()
-        self._audiothread.join()
-        self._audiothread = None
         if self._video:
             self._videothread.join()
             self._videothread = None
         self._aframectr <<= 0
         self._vframectr <<= 0
         self._video_started <<= False
+        self._rtaudio.stopStream()
 
     def isRunning(self):
         return self._running()

lasp/lasp_common.py

@@ -2,14 +2,43 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 from .wrappers import Window as wWindow
+import appdirs, os, shelve
 """
 Common definitions used throughout the code.
 """
 
-__all__ = ['P_REF', 'FreqWeighting', 'TimeWeighting', 'getTime', 'calfile',
-           'getFreq',
+__all__ = ['P_REF', 'FreqWeighting', 'TimeWeighting', 'getTime', 
+           'getFreq', 'lasp_shelve',
            'W_REF', 'U_REF', 'I_REF']
 
+lasp_appdir = appdirs.user_data_dir('Lasp', 'ASCEE')
+
+
+if not os.path.exists(lasp_appdir):
+    try:
+        os.mkdir(lasp_appdir)
+    except:
+        print('Fatal error: could not create application directory')
+        exit(1)
+
+
+class lasp_shelve:
+    refcount = 0
+    shelve = None
+
+    def __enter__(self):
+        if lasp_shelve.shelve is None:
+            assert lasp_shelve.refcount == 0
+            lasp_shelve.shelve = shelve.open(os.path.join(lasp_appdir, 'config.shelve'))
+        lasp_shelve.refcount += 1
+        return lasp_shelve.shelve
+
+    def __exit__(self, type, value, traceback):
+        lasp_shelve.refcount -= 1
+        if lasp_shelve.refcount == 0:
+            lasp_shelve.shelve.close()
+            lasp_shelve.shelve = None
+
 # Reference sound pressure level
 P_REF = 2e-5
 
@@ -19,9 +48,6 @@ I_REF = 1e-12  # 1 picoWatt/ m^2
 # Reference velocity for sound velocity level
 U_REF = 5e-8
 
-# Todo: fix This
-calfile = None
-
 
 class Window:
     hann = (wWindow.hann, 'Hann')
@@ -49,7 +75,6 @@ class Window:
     def getCurrent(cb):
         return Window.types[cb.currentIndex()]
 
-
 class TimeWeighting:
     none = (None, 'Raw (no time weighting)')
     uufast = (1e-4, '0.1 ms')
@@ -77,7 +102,6 @@ class TimeWeighting:
     def getCurrent(cb):
         return TimeWeighting.types[cb.currentIndex()]
 
-
 class FreqWeighting:
     """
     Frequency weighting types
@@ -104,7 +128,6 @@ class FreqWeighting:
     def getCurrent(cb):
         return FreqWeighting.types[cb.currentIndex()]
 
-
 def getTime(fs, N, start=0):
     """
     Return a time array for given number of points and sampling frequency.
@@ -117,7 +140,6 @@ def getTime(fs, N, start=0):
     assert N > 0 and fs > 0
     return np.linspace(start, start + N/fs, N, endpoint=False)
 
-
 def getFreq(fs, nfft):
     """
     return an array of frequencies for single-sided spectra

lasp/lasp_measurement.py

@@ -336,7 +336,7 @@ class Measurement:
             f.attrs['sensitivity'] = sens
         self._sens = sens
 
-    def exportAsWave(self, fn=None, force=False, sampwidth=None):
+    def exportAsWave(self, fn=None, force=False, newsampwidth=2, normalize=True):
         """
         Export measurement file as wave. In case the measurement data is stored
         as floats, the values are scaled between 0 and 1
@@ -348,9 +348,12 @@ class Measurement:
             force: If True, overwrites any existing files with the given name
             , otherwise a RuntimeError is raised.
 
-            sampwidth: sample width in bytes with which to export the data.
+            newsampwidth: sample width in bytes with which to export the data.
             This should only be given in case the measurement data is stored as
-            floating point values, otherwise an
+            floating point values, otherwise an error is thrown
+
+            normalize: If set: normalize the level to something sensible.
+
 
         """
         if fn is None:
@@ -362,39 +365,40 @@ class Measurement:
 
         if os.path.exists(fn) and not force:
             raise RuntimeError(f'File already exists: {fn}')
-        with self.file() as f:
-            audio = f['audio'][:]
 
-        if isinstance(audio.dtype, float):
-            if sampwidth is None:
-                raise ValueError('sampwidth parameter should be given '
-                                 'for float data in measurement file')
-            elif sampwidth == 2:
-                itype = np.int16
-            elif sampwidth == 4:
-                itype = np.int32
+        with self.file() as f:
+
+            audio = f['audio']
+            oldsampwidth = getSampWidth(audio.dtype)
+
+            max_ = 1.
+            if normalize:
+                # Find maximum value
+                for block in self.iterBlocks(f):
+                    blockmax = np.max(np.abs(block))
+                    max_ = blockmax if blockmax > max_ else max_
+                # Scale with maximum value only if we have a nonzero maximum value.
+                if max_ == 0.:
+                    max_ = 1.
+
+            if newsampwidth == 2:
+                newtype = np.int16
+            elif newsampwidth == 4:
+                newtype = np.int32
             else:
                 raise ValueError('Invalid sample width, should be 2 or 4')
 
-            # Find maximum
-            max = 0.
-            for block in self.iterBlocks():
-                blockmax = np.max(np.abs(block))
-                if blockmax > max:
-                    max = blockmax
-            # Scale with maximum value only if we have a nonzero maximum value.
-            if max == 0.:
-                max = 1.
-            scalefac = 2**(8 * sampwidth) / max
+            scalefac = 2**(8*(newsampwidth-oldsampwidth))
+            if normalize or isinstance(audio.dtype, float):
+                scalefac *= .01*max
 
-        with wave.open(fn, 'w') as wf:
-            wf.setparams((self.nchannels, self.sampwidth, self.samplerate, 0,
-                          'NONE', 'NONE'))
-            for block in self.iterBlocks():
-                if isinstance(block.dtype, float):
+            with wave.open(fn, 'w') as wf:
+                wf.setparams((self.nchannels, self.sampwidth,
+                              self.samplerate, 0, 'NONE', 'NONE'))
+                for block in self.iterBlocks(f):
                     # Convert block to integral data type
-                    block = (block * scalefac).astype(itype)
-                wf.writeframes(np.asfortranarray(block).tobytes())
+                    block = (block*scalefac).astype(newtype)
+                    wf.writeframes(np.asfortranarray(block).tobytes())
 
     @staticmethod
     def fromtxt(fn,
@@ -431,6 +435,8 @@ class Measurement:
             timestamp = os.path.getmtime(fn)
         if mfn is None:
             mfn = os.path.splitext(fn)[0] + '.h5'
+        else:
+            mfn = os.path.splitext(mfn)[0] + '.h5'
 
         dat = np.loadtxt(fn, skiprows=skiprows, delimiter=delimiter)
         if firstcoltime:

lasp/lasp_octavefilter.py

@@ -3,22 +3,68 @@
 """!
 Author: J.A. de Jong - ASCEE
 
-Provides the FIR implementation of the octave filter bank
+Provides the implementations of (fractional) octave filter banks
 
 """
-__all__ = ['OctaveFilterBank', 'ThirdOctaveFilterBank']
+__all__ = ['FirOctaveFilterBank', 'FirThirdOctaveFilterBank',
+           'OverallFilterBank', 'SosOctaveFilterBank',
+           'SosThirdOctaveFilterBank']
 
-from .filter.bandpass_fir import OctaveBankDesigner, ThirdOctaveBankDesigner
-from .wrappers import Decimator, FilterBank as pyxFilterBank
+from .filter.filterbank_design import (OctaveBankDesigner,
+                                       ThirdOctaveBankDesigner)
+from .wrappers import (Decimator, FilterBank as pyxFilterBank,
+                       SosFilterBank as pyxSosFilterBank)
 import numpy as np
 
 
-class FilterBank:
+class OverallFilterBank:
     """
-    Single channel octave filter bank implementation
+    Dummy type filter bank. Does nothing special, only returns output in a
+    sensible way
     """
 
     def __init__(self, fs):
+        """
+        Initialize overall filter bank
+        """
+        self.fs = fs
+        self.N = 0
+        self.xs = [0]
+
+    def filter_(self, data):
+        """
+        Filter input data
+        """
+        assert data.ndim == 2
+        assert data.shape[1] == 1, "invalid number of channels, should be 1"
+
+        if data.shape[0] == 0:
+            return {}
+
+        # Output given as a dictionary with x as the key
+        output = {}
+
+        tstart = self.N / self.fs
+        Ncur = data.shape[0]
+        tend = tstart + Ncur / self.fs
+
+        t = np.linspace(tstart, tend, Ncur, endpoint=False)
+        self.N += Ncur
+
+        output['Overall'] = {'t': t, 'data': data, 'x': 0}
+        return output
+
+    def decimation(self, x):
+        return [1]
+
+
+class FirFilterBank:
+    """
+    Single channel (fractional) octave filter bank implementation, based on FIR
+    filters and sample rate decimation.
+    """
+
+    def __init__(self, fs, xmin, xmax):
         """
         Initialize a OctaveFilterBank object.
 
@@ -29,7 +75,10 @@ class FilterBank:
         assert np.isclose(fs, 48000), "Only sampling frequency" \
             " available is 48 kHz"
 
-        maxdecimation = self.decimation(self.xs[0])
+        self.fs  = fs
+        self.xs = list(range(xmin, xmax + 1))
+
+        maxdecimation = self.designer.firDecimation(self.xs[0])
         self.decimators = []
         for dec in maxdecimation:
             self.decimators.append(Decimator(1, dec))
@@ -43,7 +92,7 @@ class FilterBank:
         self.filterbanks = []
         # Sort the x values in categories according to the required decimation
         for x in self.xs:
-            dec = self.decimation(x)
+            dec = self.designer.firDecimation(x)
             if len(dec) == 1 and dec[0] == 1:
                 xs_d1.append(x)
             elif len(dec) == 1 and dec[0] == 4:
@@ -60,24 +109,25 @@ class FilterBank:
         xs_all = [xs_d1, xs_d4, xs_d16, xs_d64, xs_d256]
         for xs in xs_all:
             nominals_txt = []
-            firs = np.empty((self.L, len(xs)), order='F')
-            for i, x in enumerate(xs):
-                #  These are the filters that do not require lasp_decimation
-                #  prior to filtering
-                nominals_txt.append(self.nominal_txt(x))
-                firs[:, i] = self.createFilter(fs, x)
-            filterbank = {'fb': pyxFilterBank(firs, 1024),
-                          'xs': xs,
-                          'nominals': nominals_txt}
-            self.filterbanks.append(filterbank)
+            if len(xs) > 0:
+                firs = np.empty((self.designer.firFilterLength, len(xs)), order='F')
+                for i, x in enumerate(xs):
+                    #  These are the filters that do not require lasp_decimation
+                    #  prior to filtering
+                    nominals_txt.append(self.designer.nominal_txt(x))
+                    firs[:, i] = self.designer.createFirFilter(fs, x)
+                filterbank = {'fb': pyxFilterBank(firs, 1024),
+                              'xs': xs,
+                              'nominals': nominals_txt}
+                self.filterbanks.append(filterbank)
 
         # Sample input counter.
         self.N = 0
 
-        # Filter output counters
         self.dec = [1, 4, 16, 64, 256]
 
-        # These intial delays are found experimentally using a toneburst
+        # Filter output counters
+        # These intial delays are found 'experimentally' using a toneburst
         # response.
         self.Nf = [915, 806, 780, 582, 338]
 
@@ -105,8 +155,9 @@ class FilterBank:
 
             t = np.linspace(tstart, tend, Nf, endpoint=False)
             self.Nf[dec_stage] += Nf
-            for i, nom in enumerate(self.filterbanks[dec_stage]['nominals']):
-                output[nom] = {'t': t, 'data': filtered[:, [i]]}
+            for i, nom_txt in enumerate(self.filterbanks[dec_stage]['nominals']):
+                x = self.designer.nominal_txt_tox(nom_txt)
+                output[nom_txt] = {'t': t, 'data': filtered[:, [i]], 'x': x}
 
         return output
 
@@ -122,36 +173,141 @@ class FilterBank:
 
         # Output given as a dictionary with x as the key
         output = {}
+        output_unsorted = {}
         self.N += data.shape[0]
 
-        output = {**output, **self.filterd(0, data)}
+        output_unsorted = {**output_unsorted, **self.filterd(0, data)}
 
         for i in range(len(self.decimators)):
             dec_stage = i+1
             if data.shape[0] > 0:
                 # Apply a decimation stage
                 data = self.decimators[i].decimate(data)
-                output = {**output, **self.filterd(dec_stage, data)}
+                output_unsorted = {**output_unsorted,
+                                   **self.filterd(dec_stage, data)}
+
+        # Create sorted output
+        for x in self.xs:
+            nom_txt = self.designer.nominal_txt(x)
+            output[nom_txt] = output_unsorted[nom_txt]
 
         return output
 
+    def decimation(self, x):
+        return self.designer.firDecimation(x)
 
-class OctaveFilterBank(FilterBank, OctaveBankDesigner):
+
+class FirOctaveFilterBank(FirFilterBank):
     """
     Filter bank which uses FIR filtering for each octave frequency band
     """
 
-    def __init__(self, fs):
-        OctaveBankDesigner.__init__(self)
-        FilterBank.__init__(self, fs)
+    def __init__(self, fs, xmin, xmax):
+        self.designer = OctaveBankDesigner()
+        FirFilterBank.__init__(self, fs, xmin, xmax)
 
 
-class ThirdOctaveFilterBank(FilterBank, ThirdOctaveBankDesigner):
+class FirThirdOctaveFilterBank(FirFilterBank):
     """
     Filter bank which uses FIR filtering for each one-third octave frequency
     band.
     """
 
-    def __init__(self, fs):
-        ThirdOctaveBankDesigner.__init__(self)
-        FilterBank.__init__(self, fs)
+    def __init__(self, fs, xmin, xmax):
+        self.designer = ThirdOctaveBankDesigner(fs)
+        FirFilterBank.__init__(self, fs, xmin, xmax)
+
+
+class SosFilterBank:
+    def __init__(self, fs, xmin, xmax):
+        """
+        Initialize a second order sections filterbank
+
+        Args:
+            fs: Sampling frequency [Hz]
+            xmin: Minimum value for the bands
+            xmax: Maximum value for the bands
+        """
+        self.fs = fs
+        self.xs = list(range(xmin, xmax + 1))
+        nfilt = len(self.xs)
+        self._fb = pyxSosFilterBank(nfilt, 5)
+        for i, x in enumerate(self.xs):
+            sos = self.designer.createSOSFilter(x)
+            self._fb.setFilter(i, sos)
+
+        self.xmin = xmin
+        self.xmax = xmax
+        self.N = 0
+
+    def filter_(self, data):
+        """
+        Filter input data
+        """
+        assert data.ndim == 2
+        assert data.shape[1] == 1, "invalid number of channels, should be 1"
+
+        if data.shape[0] == 0:
+            return {}
+
+        filtered_data = self._fb.filter_(data)
+
+        # Output given as a dictionary with nom_txt as the key
+        output = {}
+
+        tstart = self.N / self.fs
+        Ncur = data.shape[0]
+        tend = tstart + Ncur / self.fs
+
+        t = np.linspace(tstart, tend, Ncur, endpoint=False)
+        self.N += Ncur
+
+        for i, x in enumerate(self.xs):
+            # '31.5' to '16k'
+            nom_txt = self.designer.nominal_txt(x)
+            output[nom_txt] = {'t': t, 'data': filtered_data[:, [i]], 'x': x}
+
+        return output
+
+    def decimation(self, x):
+        return [1]
+
+
+class SosThirdOctaveFilterBank(SosFilterBank):
+    """
+    Filter bank which uses FIR filtering for each one-third octave frequency
+    band.
+    """
+
+    def __init__(self, fs, xmin, xmax):
+        """
+        Initialize a second order sections filterbank.
+
+        Args:
+            fs: Sampling frequency [Hz]
+            xmin: Minimum value for the bands
+            xmax: Maximum value for the bands
+        """
+        self.designer = ThirdOctaveBankDesigner(fs)
+        SosFilterBank.__init__(self, fs, xmin, xmax)
+
+
+class SosOctaveFilterBank(SosFilterBank):
+    """
+    Filter bank which uses FIR filtering for each one-third octave frequency
+    band.
+    """
+
+    def __init__(self, fs, xmin, xmax):
+        """
+        Initialize a second order sections filterbank.
+
+        Args:
+            fs: Sampling frequency [Hz]
+            xmin: Minimum value for the bands
+            xmax: Maximum value for the bands
+        """
+        self.designer = OctaveBankDesigner(fs)
+        SosFilterBank.__init__(self, fs, xmin, xmax)
+
+

lasp/lasp_record.py

@@ -1,8 +1,6 @@
-#!/usr/bin/env python3
+#!/usr/bin/python3
 # -*- coding: utf-8 -*-
 """
-Created on Sat Mar 10 08:28:03 2018
-
 Read data from stream and record sound and video at the same time
 """
 import numpy as np
@@ -10,11 +8,19 @@ from .lasp_atomic import Atomic
 from threading import Condition
 from .lasp_avstream import AvType, AvStream
 import h5py
+import dataclasses
+import os
 import time
 
+@dataclasses.dataclass
+class RecordStatus:
+    curT: float
+    done: bool
 
 class Recording:
-    def __init__(self, fn, stream, rectime=None):
+
+    def __init__(self, fn, stream, rectime=None, wait = True,
+                 progressCallback=None):
         """
 
         Args:
@@ -25,6 +31,7 @@ class Recording:
         ext = '.h5'
         if ext not in fn:
             fn += ext
+
         self._stream = stream
         self.blocksize = stream.blocksize
         self.samplerate = stream.samplerate
@@ -34,98 +41,130 @@ class Recording:
         self._fn = fn
 
         self._video_frame_positions = []
+        self._curT_rounded_to_seconds = 0
 
         self._aframeno = Atomic(0)
         self._vframeno = 0
 
-    def start(self):
-        stream = self._stream
+        self._progressCallback = progressCallback 
+        self._wait = wait
 
-        with h5py.File(self._fn, 'w') as f:
-            self._ad = f.create_dataset('audio',
-                                        (1, stream.blocksize, stream.nchannels),
-                                        dtype=stream.dtype,
-                                        maxshape=(None, stream.blocksize,
-                                                  stream.nchannels),
+        self._f = h5py.File(self._fn, 'w')
+        self._deleteFile = False
+
+    def setDelete(self, val: bool):
+        self._deleteFile = val
+
+    def __enter__(self):
+        """
+
+        with self._recording(wait=False):
+            event_loop_here()
+
+        or:
+
+        with Recording(wait=True):
+            pass
+        """
+
+        stream = self._stream
+        f = self._f
+
+        self._ad = f.create_dataset('audio',
+                                    (1, stream.blocksize, stream.nchannels),
+                                    dtype=stream.numpy_dtype,
+                                    maxshape=(None, stream.blocksize,
+                                              stream.nchannels),
+                                    compression='gzip'
+                                    )
+        if stream.hasVideo():
+            video_x, video_y = stream.video_x, stream.video_y
+            self._vd = f.create_dataset('video',
+                                        (1, video_y, video_x, 3),
+                                    dtype='uint8',
+                                    maxshape=(
+                                            None, video_y, video_x, 3),
                                         compression='gzip'
                                         )
-            if stream.hasVideo():
-                video_x, video_y = stream.video_x, stream.video_y
-                self._vd = f.create_dataset('video',
-                                            (1, video_y, video_x, 3),
-                                            dtype='uint8',
-                                            maxshape=(
-                                                None, video_y, video_x, 3),
-                                            compression='gzip'
-                                            )
 
-            f.attrs['samplerate'] = stream.samplerate
-            f.attrs['nchannels'] = stream.nchannels
-            f.attrs['blocksize'] = stream.blocksize
-            f.attrs['sensitivity'] = stream.sensitivity
-            f.attrs['time'] = time.time()
-            self._running <<= True
-            # Videothread is going to start
+        f.attrs['samplerate'] = stream.samplerate
+        f.attrs['nchannels'] = stream.nchannels
+        f.attrs['blocksize'] = stream.blocksize
+        f.attrs['sensitivity'] = stream.sensitivity
+        f.attrs['time'] = time.time()
+        self._running <<= True
 
-            if not stream.isRunning():
-                stream.start()
+        if not stream.isRunning():
+            stream.start()
 
-            stream.addCallback(self._callback)
+        print('Starting record....')
+        stream.addCallback(self._aCallback, AvType.audio_input)
+        if stream.hasVideo():
+            stream.addCallback(self._aCallback, AvType.audio_input)
+
+        if self._wait:
             with self._running_cond:
+                print('Stop recording with CTRL-C')
                 try:
-                    print('Starting record....')
                     while self._running:
                         self._running_cond.wait()
                 except KeyboardInterrupt:
                     print("Keyboard interrupt on record")
                     self._running <<= False
 
-            stream.removeCallback(self._callback)
 
-            if stream.hasVideo():
-                f['video_frame_positions'] = self._video_frame_positions
-
-            print('\nEnding record')
-
-    def stop(self):
+    def __exit__(self, type, value, traceback):
         self._running <<= False
-        with self._running_cond:
-            self._running_cond.notify()
+        stream = self._stream
+        stream.removeCallback(self._aCallback, AvType.audio_input)
+        if stream.hasVideo():
+            stream.removeCallback(self._vCallback, AvType.video_input)
+            f['video_frame_positions'] = self._video_frame_positions
 
-    def _callback(self, _type, data, aframe, vframe):
-        if not self._stream.isRunning():
-            self._running <<= False
-            with self._running_cond:
-                self._running_cond.notify()
+        self._f.close()
+        print('\nEnding record')
+        if self._deleteFile:
+            try:
+                os.remove(self._fn)
+            except Exception as e:
+                print(f'Error deleting file: {self._fn}')
 
-        if _type == AvType.audio:
-            self._aCallback(data, aframe)
-        elif _type == AvType.video:
-            self._vCallback(data)
 
-    def _aCallback(self, frames, aframe):
-        print('.', end='')
+
+    def _aCallback(self, indata, outdata, aframe):
+
         curT = self._aframeno()*self.blocksize/self.samplerate
+        recstatus = RecordStatus(
+            curT = curT,
+            done = False)
+        if self._progressCallback is not None:
+            self._progressCallback(recstatus)
+
+        curT_rounded_to_seconds = int(curT)
+        if curT_rounded_to_seconds > self._curT_rounded_to_seconds:
+            self._curT_rounded_to_seconds = curT_rounded_to_seconds
+            print(f'{curT_rounded_to_seconds}', end='', flush=True)
+        else:
+            print('.', end='', flush=True)
+
         if self.rectime is not None and curT > self.rectime:
             # We are done!
             self._running <<= False
             with self._running_cond:
                 self._running_cond.notify()
-                return
+            if self._progressCallback is not None:
+                recstatus.done = True
+                self._progressCallback(recstatus)
+            return
 
         self._ad.resize(self._aframeno()+1, axis=0)
-        self._ad[self._aframeno(), :, :] = frames
+        self._ad[self._aframeno(), :, :] = indata
         self._aframeno += 1
 
-    def _vCallback(self, frame):
+    def _vCallback(self, frame, framectr):
         self._video_frame_positions.append(self._aframeno())
         vframeno = self._vframeno
         self._vd.resize(vframeno+1, axis=0)
         self._vd[vframeno, :, :] = frame
         self._vframeno += 1
 
-
-if __name__ == '__main__':
-    stream = AvStream()
-    rec = Recording('test', stream, 5)
-    rec.start()

lasp/lasp_slm.py

@@ -47,7 +47,7 @@ class SLM:
         self._Lmax = 0.
 
         # Storage for computing the equivalent level
-        self._sq = 0.
+        self._sq = 0. # Square of the level data, storage
         self._N = 0
         self._Leq = 0.
 

lasp/lasp_weighcal.py

@@ -4,10 +4,8 @@
 Weighting and calibration filter in one
 @author: J.A. de Jong - ASCEE
 """
-from .filter.freqweighting_fir import A, C
 from .lasp_common import FreqWeighting
-from .filter.fir_design import (arbitrary_fir_design,
-                                freqResponse as frp)
+from .filter import (A, C, arbitrary_fir_design, freqResponse as frp)
 from lasp.lasp_config import ones, empty
 from .wrappers import FilterBank
 import numpy as np

lasp/tools/plot.py

@@ -5,7 +5,7 @@ Author: J.A. de Jong - ASCEE
 
 Description:
 """
-__all__ = ['close', 'Figure', 'Bode', 'PS', 'PSD']
+__all__ = ['Figure', 'Bode', 'PS', 'PSD']
 
 from .config import getReportQuality
 import matplotlib.pyplot as plt

lasp/wrappers.pyx

@@ -2,8 +2,70 @@
 This file contains the Cython wrapper functions to C implementations.
 """
 include "config.pxi"
+IF LASP_FLOAT == "double":
+    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
 
-__all__ = ['AvPowerSpectra']
+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']
 
 setTracerLevel(15)
 cdef extern from "cblas.h":
@@ -65,11 +127,11 @@ cdef class Fft:
                                         timedata.shape[1],
                                         &timedata[0,0],
                                         False)
+        with nogil:
+            Fft_fft(self._fft,&t,&r)
 
-        Fft_fft(self._fft,&t,&r)
-
-        dmat_free(&t)
-        cmat_free(&r)
+            dmat_free(&t)
+            cmat_free(&r)
 
         return result
 
@@ -97,10 +159,11 @@ cdef class Fft:
                                         &timedata_view[0,0],
                                         False)
 
-        Fft_ifft(self._fft,&f,&t)
+        with nogil:
+            Fft_ifft(self._fft,&f,&t)
 
-        dmat_free(&t)
-        cmat_free(&f)
+            dmat_free(&t)
+            cmat_free(&f)
 
         return timedata
 
@@ -127,7 +190,7 @@ cdef extern from "lasp_ps.h":
 
     void PowerSpectra_compute(const c_PowerSpectra* ps,
                              const dmat * timedata,
-                             cmat * result)
+                             cmat * result) nogil
 
 
     void PowerSpectra_free(c_PowerSpectra*)
@@ -174,8 +237,8 @@ cdef class PowerSpectra:
                                        False)
 
 
-
-        PowerSpectra_compute(self._ps,&td,&result_mat)
+        with nogil:
+            PowerSpectra_compute(self._ps,&td,&result_mat)
 
         dmat_free(&td)
         cmat_free(&result_mat)
@@ -197,7 +260,7 @@ cdef extern from "lasp_aps.h":
                                            const vd* weighting)
 
     cmat* AvPowerSpectra_addTimeData(const c_AvPowerSpectra* ps,
-                                     const dmat * timedata)
+                                     const dmat * timedata) nogil
 
 
     void AvPowerSpectra_free(c_AvPowerSpectra*)
@@ -257,16 +320,6 @@ cdef class AvPowerSpectra:
                                &timedata[0,0],
                                False)
 
-        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]
-
         result = np.empty((self.nfft//2+1,nchannels,nchannels),
                           dtype = NUMPY_COMPLEX_TYPE,
                           order='F')
@@ -277,45 +330,111 @@ cdef class AvPowerSpectra:
                                           nchannels*nchannels,
                                           &result_view[0,0,0],
                                           False)
-        # Copy result
-        cmat_copy(&res,result_ptr)
+        with nogil:
+            result_ptr = AvPowerSpectra_addTimeData(self.aps,
+                                                &td)
 
-        cmat_free(&res)
-        dmat_free(&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_filterbank.h":
-    ctypedef struct c_FilterBank "FilterBank"
-    c_FilterBank* FilterBank_create(const dmat* h,const us nfft) nogil
-    dmat FilterBank_filter(c_FilterBank* fb,const vd* x) nogil
-    void FilterBank_free(c_FilterBank* fb) nogil
+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_FilterBank* fb
+        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 = FilterBank_create(&hmat,nfft)
+        self.fb = Firfilterbank_create(&hmat,nfft)
         dmat_free(&hmat)
         if not self.fb:
             raise RuntimeError('Error creating FilberBank')
 
     def __dealloc__(self):
         if self.fb:
-            FilterBank_free(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 = FilterBank_filter(self.fb,&input_vd)
+        
+        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 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
+    def __cinit__(self,const us filterbank_size, const us nsections):
+
+        self.fb = Sosfilterbank_create(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
+
+        """
+        cdef dmat coefs = dmat_foreign_data(sos.size,1,
+                                            &sos[0, 0],False)
+        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)
 
         # Steal the pointer from output
         result = dmat_to_ndarray(&output,True)
@@ -400,3 +519,54 @@ cdef class SPLowpass:
         vd_free(&input_vd)
 
         return result
+
+
+cdef extern from "lasp_siggen.h":
+    ctypedef struct c_Siggen "Siggen"
+    c_Siggen* Siggen_Whitenoise_create(d fs, d level_dB)
+    c_Siggen* Siggen_Sinewave_create(d fs, d freq, d level_dB)
+    void Siggen_genSignal(c_Siggen*, vd* samples) nogil
+    void Siggen_free(c_Siggen*)
+
+
+cdef class Siggen:
+    cdef c_Siggen *_siggen
+
+    def __cinit__(self):
+        self._siggen = NULL
+
+    def __dealloc__(self):
+        if self._siggen:
+            Siggen_free(self._siggen)
+
+    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
+
+
+    @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 whiteNoise(d fs, d level_dB):
+        cdef c_Siggen* c_siggen = Siggen_Whitenoise_create(fs, level_dB)
+        siggen = Siggen()
+        siggen._siggen = c_siggen
+        return siggen

scripts/lasp_calibrate.py

@@ -26,7 +26,7 @@ parser.add_argument('--gain-setting', '-g',
                     type=float, default=gain_default)
 
 parser.add_argument(
-    'fn', help='File name of calibration measurement', type=str)
+    'fn', help='File name of calibration measurement', type=str, default=None)
 
 parser.add_argument('--channel', help='Channel of the device to calibrate, '
                     + 'default = 0',
@@ -51,10 +51,11 @@ prms = P_REF*10**(args.spl/20)
 sens = Vrms / prms
 
 print(f'Computed sensitivity: {sens[args.channel]:.5} V/Pa')
-print('Searching for files in directory to apply sensitivity value to...')
-dir_ = os.path.dirname(args.fn)
-for f in os.listdir(dir_):
-    yn = input(f'Apply sensitivity to {f}? [Y/n]')
-    if yn in ['', 'Y', 'y']:
-        meas = Measurement(os.path.join(dir_, f))
-        meas.sensitivity = sens
+if args.fn:
+    print('Searching for files in directory to apply sensitivity value to...')
+    dir_ = os.path.dirname(args.fn)
+    for f in os.listdir(dir_):
+        yn = input(f'Apply sensitivity to {f}? [Y/n]')
+        if yn in ['', 'Y', 'y']:
+            meas = Measurement(os.path.join(dir_, f))
+            meas.sensitivity = sens

scripts/lasp_record

@@ -1,5 +1,8 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 import argparse
+
+
+
 parser = argparse.ArgumentParser(
     description='Acquire data and store a measurement file'
 )
@@ -8,38 +11,49 @@ parser.add_argument('filename', type=str,
                     ' Extension is automatically added.')
 parser.add_argument('--duration', '-d', type=float,
                     help='The recording duration in [s]')
-parser.add_argument('--comment', '-c', type=str,
-                    help='Add a measurement comment, optionally')
 
 device_help = 'DAQ Device to record from'
 parser.add_argument('--input-daq', '-i', help=device_help, type=str,
-                    choices=['roga', 'umik', 'nidaq', 'default'], default='roga')
+                    default='Default')
 
 args = parser.parse_args()
 
-device_str = args.input_daq
-
-if device_str == 'nidaq':
-   
-    # Not-integrated way to record with the NI USB 4431 DAQ device
-    from lasp.device.record_ni import USBDAQRecording
-    rec =  USBDAQRecording(args.filename, [0, 1])
-    rec.start()
-    exit(0)
-
 
+from lasp.lasp_avstream import AvStream, AvType
 from lasp.lasp_record import Recording
-from lasp.lasp_avstream import AvStream
-from lasp.device.lasp_daqconfig import default_soundcard, roga_plugndaq, umik
+from lasp.device import DAQConfiguration, RtAudio
 
-if 'roga' == device_str:
-    device = roga_plugndaq
-elif 'default' == device_str:
-    device = default_soundcard
-elif 'umik' == device_str:
-    device = umik
+config = DAQConfiguration.loadConfigs()[args.input_daq]
+
+print(config)
+rtaudio = RtAudio()
+count = rtaudio.getDeviceCount()
+devices = [rtaudio.getDeviceInfo(i) for i in range(count)]
+
+input_devices = {}
+for device in devices:
+    if device.inputchannels >= 0:
+        input_devices[device.name] = device
+
+try:
+    input_device = input_devices[config.input_device_name]
+except KeyError:
+    raise RuntimeError(f'Input device {config.input_device_name} not available')
+
+print(input_device)
+
+stream = AvStream(input_device,
+                  AvType.audio_input,
+                  config)
 
-stream = AvStream(device)
 rec = Recording(args.filename, stream, args.duration)
-rec.start()
+with rec:
+    pass
+
+print('Stopping stream...')
 stream.stop()
+
+print('Stream stopped')
+print('Closing stream...')
+stream.close()
+print('Stream closed')

scripts/lasp_siggen

@@ -0,0 +1,66 @@
+#!/usr/bin/python3
+import argparse
+import numpy as np
+
+
+parser = argparse.ArgumentParser(
+    description='Play a sine wave'
+)
+device_help = 'DAQ Device to play to'
+parser.add_argument('--device', '-d', help=device_help, type=str,
+                    default='Default')
+
+args = parser.parse_args()
+
+from lasp.lasp_avstream import AvStream, AvType
+from lasp.device import DAQConfiguration, RtAudio
+
+config = DAQConfiguration.loadConfigs()[args.device]
+
+rtaudio = RtAudio()
+count = rtaudio.getDeviceCount()
+devices = [rtaudio.getDeviceInfo(i) for i in range(count)]
+
+output_devices = {}
+for device in devices:
+    if device.outputchannels >= 0:
+        output_devices[device.name] = device
+
+try:
+    output_device = output_devices[config.output_device_name]
+except KeyError:
+    raise RuntimeError(f'output device {config.output_device_name} not available')
+
+samplerate = int(config.en_output_rate)
+stream = AvStream(output_device,
+                  AvType.audio_output,
+                  config)
+
+# freq = 440.
+freq = 1000.
+omg = 2*np.pi*freq
+
+
+def mycallback(indata, outdata, blockctr):
+    frames = outdata.shape[0]
+    nchannels = outdata.shape[1]
+    # nchannels = 1
+    streamtime = blockctr*frames/samplerate
+    t = np.linspace(streamtime, streamtime + frames/samplerate,
+            frames)[np.newaxis, :]
+    outp = 0.01*np.sin(omg*t)
+    for i in range(nchannels):
+        outdata[:,i] = ((2**16-1)*outp).astype(np.int16)
+
+stream.addCallback(mycallback, AvType.audio_output)
+stream.start()
+
+input()
+
+print('Stopping stream...')
+stream.stop()
+
+print('Stream stopped')
+print('Closing stream...')
+stream.close()
+print('Stream closed')

test/fft_test.py

@@ -6,7 +6,7 @@ Created on Mon Jan 15 19:45:33 2018
 @author: anne
 """
 import numpy as np
-from lasp import Fft
+from lasp.wrappers import Fft
 
 nfft=9
 print('nfft:',nfft)

test/filterbank_test.py

@@ -33,4 +33,4 @@ print('res shape %i:\n' %res.shape[0])
 print(res)
 #for i in range(4):
 #    res = fb.filter_(data)
-#    print('res:\n',res)
+#    print('res:\n',res)

test/sosfilterbank_test.py

@@ -0,0 +1,26 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Jan 15 19:45:33 2018
+
+@author: anne
+"""
+import numpy as np
+from lasp.wrappers import SosFilterBank
+
+fb = SosFilterBank(1, 2)
+
+newfilter1 = np.array([0,1.,0.,1.,0.,0.])[:,np.newaxis] # Unit sample delay
+newfilter2 = np.array([0,1.,0.,1.,0.,0.])[:,np.newaxis] # Unit sample delay
+newfilter = np.vstack((newfilter1, newfilter2))
+
+fb.setFilter(0, newfilter)
+x = np.zeros(10)
+x[5]=1
+x[8] = 3
+x = x[:,np.newaxis]
+y = fb.filter_(x)
+print(x)
+print(y)
+y = fb.filter_(x)
+print(y)

test/test_rtaudio.py

@@ -0,0 +1,35 @@
+#!/usr/bin/python3
+import numpy as np
+from lasp_rtaudio import RtAudio, SampleFormat, Format_SINT32, Format_FLOAT64
+import time
+
+
+nframes = 0
+samplerate = 48000
+omg = 2*np.pi*1000
+
+def mycallback(input_, nframes, streamtime):
+    t = np.linspace(streamtime, streamtime + nframes/samplerate,
+            nframes)[np.newaxis,:]
+    outp = 0.1*np.sin(omg*t)
+    return outp, 0
+
+if __name__ == '__main__':
+    pa = RtAudio()
+    count = pa.getDeviceCount()
+    # dev = pa.getDeviceInfo(0)
+    for i in range(count):
+        dev = pa.getDeviceInfo(i)
+        print(dev)
+   
+    outputparams = {'deviceid': 0, 'nchannels': 1, 'firstchannel': 0}
+    pa.openStream(outputparams, None , Format_FLOAT64,samplerate, 512, mycallback)
+    pa.startStream()
+
+    input()
+
+    pa.stopStream()
+    pa.closeStream()
+
+
+