From 10749137eccebe4a434e3a82c9668f8dfdf46508 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Sat, 3 Sep 2022 16:10:12 +0200 Subject: [PATCH] Removed old C-code. This is not the way forward anymore --- src/lasp/c/CMakeLists.txt | 29 -- src/lasp/c/dec_filter_4.c | 128 ------- src/lasp/c/lasp_alg.c | 223 ------------ src/lasp/c/lasp_alg.h | 192 ---------- src/lasp/c/lasp_alloc.h | 33 -- src/lasp/c/lasp_aps.c | 291 ---------------- src/lasp/c/lasp_aps.h | 85 ----- src/lasp/c/lasp_assert.c | 29 -- src/lasp/c/lasp_assert.h | 44 --- src/lasp/c/lasp_bem.c | 20 -- src/lasp/c/lasp_decimation.c | 221 ------------ src/lasp/c/lasp_decimation.h | 54 --- src/lasp/c/lasp_dfifo.c | 145 -------- src/lasp/c/lasp_dfifo.h | 73 ---- src/lasp/c/lasp_eq.c | 70 ---- src/lasp/c/lasp_eq.h | 59 ---- src/lasp/c/lasp_fft.h | 109 ------ src/lasp/c/lasp_firfilterbank.c | 174 ---------- src/lasp/c/lasp_firfilterbank.h | 59 ---- src/lasp/c/lasp_mat.c | 54 --- src/lasp/c/lasp_mat.h | 597 -------------------------------- src/lasp/c/lasp_math_raw.c | 152 -------- src/lasp/c/lasp_math_raw.h | 480 ------------------------- src/lasp/c/lasp_mq.c | 367 -------------------- src/lasp/c/lasp_mq.h | 78 ----- src/lasp/c/lasp_nprocs.c | 18 - src/lasp/c/lasp_nprocs.h | 25 -- src/lasp/c/lasp_ps.c | 181 ---------- src/lasp/c/lasp_ps.h | 61 ---- src/lasp/c/lasp_pyarray.h | 108 ------ src/lasp/c/lasp_python.h | 59 ---- src/lasp/c/lasp_signals.h | 32 -- src/lasp/c/lasp_slm.c | 293 ---------------- src/lasp/c/lasp_slm.h | 96 ----- src/lasp/c/lasp_sosfilterbank.c | 264 -------------- src/lasp/c/lasp_sosfilterbank.h | 81 ----- src/lasp/c/lasp_worker.c | 206 ----------- src/lasp/c/lasp_worker.h | 65 ---- 38 files changed, 5255 deletions(-) delete mode 100644 src/lasp/c/CMakeLists.txt delete mode 100644 src/lasp/c/dec_filter_4.c delete mode 100644 src/lasp/c/lasp_alg.c delete mode 100644 src/lasp/c/lasp_alg.h delete mode 100644 src/lasp/c/lasp_alloc.h delete mode 100644 src/lasp/c/lasp_aps.c delete mode 100644 src/lasp/c/lasp_aps.h delete mode 100644 src/lasp/c/lasp_assert.c delete mode 100644 src/lasp/c/lasp_assert.h delete mode 100644 src/lasp/c/lasp_bem.c delete mode 100644 src/lasp/c/lasp_decimation.c delete mode 100644 src/lasp/c/lasp_decimation.h delete mode 100644 src/lasp/c/lasp_dfifo.c delete mode 100644 src/lasp/c/lasp_dfifo.h delete mode 100644 src/lasp/c/lasp_eq.c delete mode 100644 src/lasp/c/lasp_eq.h delete mode 100644 src/lasp/c/lasp_fft.h delete mode 100644 src/lasp/c/lasp_firfilterbank.c delete mode 100644 src/lasp/c/lasp_firfilterbank.h delete mode 100644 src/lasp/c/lasp_mat.c delete mode 100644 src/lasp/c/lasp_mat.h delete mode 100644 src/lasp/c/lasp_math_raw.c delete mode 100644 src/lasp/c/lasp_math_raw.h delete mode 100644 src/lasp/c/lasp_mq.c delete mode 100644 src/lasp/c/lasp_mq.h delete mode 100644 src/lasp/c/lasp_nprocs.c delete mode 100644 src/lasp/c/lasp_nprocs.h delete mode 100644 src/lasp/c/lasp_ps.c delete mode 100644 src/lasp/c/lasp_ps.h delete mode 100644 src/lasp/c/lasp_pyarray.h delete mode 100644 src/lasp/c/lasp_python.h delete mode 100644 src/lasp/c/lasp_signals.h delete mode 100644 src/lasp/c/lasp_slm.c delete mode 100644 src/lasp/c/lasp_slm.h delete mode 100644 src/lasp/c/lasp_sosfilterbank.c delete mode 100644 src/lasp/c/lasp_sosfilterbank.h delete mode 100644 src/lasp/c/lasp_worker.c delete mode 100644 src/lasp/c/lasp_worker.h diff --git a/src/lasp/c/CMakeLists.txt b/src/lasp/c/CMakeLists.txt deleted file mode 100644 index 905866d..0000000 --- a/src/lasp/c/CMakeLists.txt +++ /dev/null @@ -1,29 +0,0 @@ -if(!LASP_DEBUG) -SET_SOURCE_FILES_PROPERTIES(lasp_sosfilterbank.c PROPERTIES COMPILE_FLAGS -O3) -SET_SOURCE_FILES_PROPERTIES(lasp_slm.c PROPERTIES COMPILE_FLAGS -O3) -SET_SOURCE_FILES_PROPERTIES(lasp_eq.c PROPERTIES COMPILE_FLAGS -O3) -endif(!LASP_DEBUG) - - -# add_library(lasp_lib -# lasp_fft.c -# lasp_mat.c -# # lasp_math_raw.c -# # lasp_alg.c -# lasp_assert.c -# lasp_tracer.c -# # lasp_window.c -# lasp_aps.c -# lasp_ps.c -# # lasp_mq.c -# lasp_worker.c -# lasp_nprocs.c -# lasp_dfifo.c -# lasp_firfilterbank.c -# lasp_sosfilterbank.c -# lasp_decimation.c -# lasp_slm.c -# lasp_eq.c -# ) - -# target_link_libraries(lasp_lib ${FFT_LIBRARIES} ${BLAS_LIBRARIES}) diff --git a/src/lasp/c/dec_filter_4.c b/src/lasp/c/dec_filter_4.c deleted file mode 100644 index 8f11524..0000000 --- a/src/lasp/c/dec_filter_4.c +++ /dev/null @@ -1,128 +0,0 @@ -{0.000000000000000000e+00, - -2.185738981122740943e-06, - -2.794005403671438005e-06, - 9.474117030816240549e-06, - 4.047683164318276243e-05, - 8.159721257084205745e-05, - 1.081272986403995632e-04, - 8.819754747279720866e-05, - 1.292089375486244570e-18, - -1.501119890529942897e-04, - -3.184064926464186956e-04, - -4.310893684346957895e-04, - -4.094689334832080671e-04, - -2.059289266182427725e-04, - 1.632175892397718904e-04, - 6.039349831684049600e-04, - 9.609848668279323634e-04, - 1.065386827404308026e-03, - 8.003236325878445483e-04, - 1.624425171981582980e-04, - -7.076130570048387320e-04, - -1.544654803296158230e-03, - -2.031572983126828796e-03, - -1.908002710205690703e-03, - -1.080539502548179872e-03, - 3.071493710103571440e-04, - 1.878468916273352439e-03, - 3.114730098982527052e-03, - 3.515061533711236544e-03, - 2.779296151548805299e-03, - 9.515559570652546090e-04, - -1.533803816637046430e-03, - -3.936604506597823384e-03, - -5.417051906137869584e-03, - -5.310758368959681876e-03, - -3.387308238052255758e-03, - -2.020189585846376504e-17, - 3.941575036075541973e-03, - 7.193441582498504364e-03, - 8.546931560120834062e-03, - 7.242589074154118407e-03, - 3.295234447391589671e-03, - -2.391580800508445910e-03, - -8.190127435210622919e-03, - -1.217754814426674249e-02, - -1.272615168819673029e-02, - -9.085393051846636994e-03, - -1.766149937181874050e-03, - 7.423829658638344230e-03, - 1.575380997299343638e-02, - 2.029368152419728719e-02, - 1.881247426920701696e-02, - 1.060259142796622991e-02, - -3.026237731635000611e-03, - -1.877030360461057201e-02, - -3.192876928283294724e-02, - -3.747119479720732033e-02, - -3.132999667548151679e-02, - -1.158792905301173765e-02, - 2.076675004965774715e-02, - 6.175096553249778686e-02, - 1.050334155991766993e-01, - 1.431963133103141828e-01, - 1.693251638559906402e-01, - 1.785441122547395953e-01, - 1.691180120289761668e-01, - 1.428459416149487626e-01, - 1.046477658307891495e-01, - 6.144841153891738433e-02, - 2.063940827275027520e-02, - -1.150252035146304662e-02, - -3.106003822342199780e-02, - -3.710127967927729503e-02, - -3.157313813741360886e-02, - -1.853722959527762462e-02, - -2.984746593388511431e-03, - 1.044334023900628412e-02, - 1.850493210242899755e-02, - 1.993456916585750749e-02, - 1.545344532020422393e-02, - 7.271929013423104535e-03, - -1.727500581384247688e-03, - -8.873383043962101632e-03, - -1.241029376178266405e-02, - -1.185679908045850044e-02, - -7.961640663049178793e-03, - -2.321032949510159950e-03, - 3.192604628017779635e-03, - 7.004731963719125487e-03, - 8.251271964876608425e-03, - 6.931579083478471744e-03, - 3.790698415906206074e-03, - -1.938928732007791559e-17, - -3.244201578800122790e-03, - -5.075193378062196545e-03, - -5.164847959311490676e-03, - -3.744259045513907962e-03, - -1.455155744093486990e-03, - 9.003467998605640694e-04, - 2.622284972952067337e-03, - 3.306540383383136175e-03, - 2.920618605100662978e-03, - 1.755414747876200702e-03, - 2.859847735621343567e-04, - -1.002153516667671955e-03, - -1.762137666026810951e-03, - -1.867706509456178859e-03, - -1.413021131661008033e-03, - -6.438025163091647607e-04, - 1.469135513812169238e-04, - 7.190488353391707912e-04, - 9.501796827231699877e-04, - 8.500266304384044257e-04, - 5.292428644636570966e-04, - 1.415165912868260117e-04, - -1.763689026852747436e-04, - -3.456944177857647653e-04, - -3.578084944368887165e-04, - -2.589144427740892114e-04, - -1.190202164875697825e-04, - 9.922942950529402484e-19, - 6.497054701652780085e-05, - 7.525727688871557033e-05, - 5.231825144807937484e-05, - 2.280076622100770170e-05, - 4.215016359691300044e-06, - -6.989289499831618075e-07, - -0.000000000000000000e+00} diff --git a/src/lasp/c/lasp_alg.c b/src/lasp/c/lasp_alg.c deleted file mode 100644 index 6947b3e..0000000 --- a/src/lasp/c/lasp_alg.c +++ /dev/null @@ -1,223 +0,0 @@ -// lasp_alg.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// (Linear) algebra routine implementations -////////////////////////////////////////////////////////////////////// -#include "lasp_alg.h" - -void cmv_dot(const cmat* A,const vc* restrict x,vc* restrict b){ - dbgassert(A && x && b,NULLPTRDEREF); - assert_vx(b); - assert_vx(x); - dbgassert(A->n_cols == x->n_rows,SIZEINEQUAL); - dbgassert(A->n_rows == b->n_rows,SIZEINEQUAL); - - #if LASP_USE_BLAS == 1 - dbgassert(false,"Untested function. Is not functional for strides"); - /* typedef enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER; */ - /* typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113, CblasConjNoTrans=114} CBLAS_TRANSPOSE; */ - /* - void cblas_zgemv(OPENBLAS_CONST enum CBLAS_ORDER order, - OPENBLAS_CONST enum CBLAS_TRANSPOSE trans, - OPENBLAS_CONST blasint m, - OPENBLAS_CONST blasint n, - OPENBLAS_CONST double *alpha, - OPENBLAS_CONST double *a, - OPENBLAS_CONST blasint lda, - OPENBLAS_CONST double *x, - OPENBLAS_CONST blasint incx, - OPENBLAS_CONST double *beta, - double *y, - OPENBLAS_CONST blasint incy); - */ - c alpha = 1.0; - c beta = 0.0; - cblas_zgemv(CblasColMajor, - CblasNoTrans, - A->n_rows, - A->n_cols, - (d*) &alpha, /* alpha */ - (d*) A->_data, /* A */ - A->n_rows, /* lda */ - (d*) x->_data, /* */ - 1, - (d*) &beta, /* beta */ - (d*) b->_data, - 1); - - - - #else - size_t i,j; - - vc_set(b,0.0); - - iVARTRACE(20,A->n_cols); - iVARTRACE(20,A->n_rows); - - for(j=0;jn_cols;j++){ - for(i=0;in_rows;i++) { - - c* Aij = getcmatval(A,i,j); - b->_data[i] += *Aij * *getvcval(x,j); - - } - - } - - - #endif -} - -/// The code below here is not yet worked on. Should be improved to -/// directly couple to openblas, instead of using lapacke.h -#if 0 - -/* These functions can be directly linked to openBLAS */ -#define lapack_complex_double double _Complex -#define lapack_complex_float float _Complex - -#define LAPACK_ROW_MAJOR 101 -#define LAPACK_COL_MAJOR 102 - -#define LAPACK_WORK_MEMORY_ERROR -1010 -#define LAPACK_TRANSPOSE_MEMORY_ERROR -1011 - -typedef int lapack_int; - -int LAPACKE_cgelss( int matrix_layout, int m, int n, - int nrhs, lapack_complex_float* a, - int lda, lapack_complex_float* b, - int ldb, float* s, float rcond, - int* rank ); -int LAPACKE_zgelss( int matrix_layout, int m, int n, - int nrhs, lapack_complex_double* a, - int lda, lapack_complex_double* b, - int ldb, double* s, double rcond, - int* rank ); - -lapack_int LAPACKE_zgels( int matrix_layout, char trans, lapack_int m, - lapack_int n, lapack_int nrhs, - lapack_complex_double* a, lapack_int lda, - lapack_complex_double* b, lapack_int ldb ); - - - - -#if LASP_FLOAT == 64 - -#define lapack_gelss LAPACKE_zgelss -#define lapack_gels LAPACKE_zgels -#else - -#define lapack_gelss LAPACKE_cgelss -#endif - -#define max(a,b) ((a)>(b)?(a):(b)) - - -/* int lsq_solve(const cmat* A,const vc* b,vc* x){ */ - -/* int rv; */ -/* /\* M: number of rows of matrix *\/ */ -/* /\* N: Number of columns *\/ */ -/* /\* Norm: L2|b-A*x| *\/ */ -/* /\* NRHS: Number of right hand sides: Number of columns of matrix B *\/ */ - -/* assert(A->n_rows>=A->n_cols); */ -/* assert(x->size == A->n_cols); */ -/* assert(b->size == A->n_rows); */ - -/* int info; */ - -/* size_t lda = max(1,A->n_rows); */ -/* size_t ldb = max(lda,A->n_cols); */ - -/* /\* Make explicit copy of matrix A data, as it will be overwritten */ -/* * by lapack_gels *\/ */ -/* c* A_data = Pool_allocatec(&lsq_solve_pool,A->n_rows*A->n_cols); */ -/* c_copy(A_data,A->data,A->n_cols*A->n_rows); */ - -/* c* work_data = Pool_allocatec(&lsq_solve_pool,b->size); */ -/* c_copy(work_data,b->data,b->size); */ - -/* /\* Lapack documentation says: *\/ */ -/* /\* if TRANS = 'N' and m >= n, rows 1 to n of B contain the least */ -/* squares solution vectors; the residual sum of squares for the */ -/* solution in each column is given by the sum of squares of the */ -/* modulus of elements N+1 to M in that column; */ -/* *\/ */ - - -/* /\* We always assume one RHS column *\/ */ -/* const int nrhs = 1; */ - -/* /\* General Least Squares Solve *\/ */ -/* info = lapack_gels(LAPACK_COL_MAJOR, /\* Column-major ordering *\/ */ -/* 'N', */ -/* A->n_rows, /\* Number of rows in matrix *\/ */ -/* A->n_cols, /\* Number of columns *\/ */ -/* nrhs, /\* nrhs, which is number_mics *\/ */ -/* A_data, /\* The A-matrix *\/ */ -/* lda, /\* lda: the leading dimension of matrix A *\/ */ -/* work_data, /\* The b-matrix *\/ */ -/* ldb); /\* ldb: the leading dimension of b: max(1,M,N) *\/ */ - -/* if(info==0){ */ -/* c_copy(x->data,work_data,x->size); */ -/* rv = SUCCESS; */ -/* } */ -/* else { */ -/* memset(x->data,0,x->size); */ -/* WARN("LAPACK INFO VALUE"); */ -/* printf("%i\n", info ); */ -/* TRACE(15,"Solving least squares problem failed\n"); */ - -/* rv = FAILURE; */ -/* } */ - -/* return rv; */ - -/* } */ - -/* d c_normdiff(const cmat* A,const cmat* B) { */ - -/* TRACE(15,"c_normdif"); */ - -/* dbgassert(A->n_cols==B->n_cols,"Number of columns of A and B " */ -/* "should be equal"); */ -/* dbgassert(A->n_rows==B->n_rows,"Number of rows of A and B " */ -/* "should be equal"); */ - -/* size_t size = A->n_cols*A->n_rows; */ - -/* vc diff_temp = vc_al[MAX_MATRIX_SIZE]; */ - -/* c_copy(diff_temp,A->data,size); */ - -/* c alpha = -1.0; */ - -/* /\* This routine computes y <- alpha*x + beta*y *\/ */ - - -/* /\* void cblas_zaxpy(OPENBLAS_CONST blasint n, *\/ */ -/* /\* OPENBLAS_CONST double *alpha, *\/ */ -/* /\* OPENBLAS_CONST double *x, *\/ */ -/* /\* OPENBLAS_CONST blasint incx, *\/ */ -/* /\* double *y, *\/ */ -/* /\* OPENBLAS_CONST blasint incy); *\/ */ - -/* cblas_zaxpy(size, */ -/* (d*) &alpha, */ -/* (d*) B->data, */ -/* 1, */ -/* (d*) diff_temp, */ -/* 1 ); */ - -/* return c_norm(diff_temp,size); */ -/* } */ -#endif /* if 0 */ - -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_alg.h b/src/lasp/c/lasp_alg.h deleted file mode 100644 index bc41074..0000000 --- a/src/lasp/c/lasp_alg.h +++ /dev/null @@ -1,192 +0,0 @@ -// lasp_alg.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// (Linear) algebra routines on matrices and vectors -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_ALG_H -#define LASP_ALG_H -#include "lasp_config.h" -#include "lasp_mat.h" - -/** - * Compute the dot product of two vectors of floats - * - * @param a First vector - * @param b Second second vector - * @return dot product as float - */ -static inline d vd_dot(const vd * a,const vd* b) { - dbgassert(a && b,NULLPTRDEREF); - assert_vx(a); - assert_vx(b); - dbgassert(a->n_rows == b->n_rows,SIZEINEQUAL); - return d_dot(getvdval(a,0),getvdval(b,0),a->n_rows); -} - -/** - * y = fac * y - * - * @param y - * @param fac scale factor - */ -static inline void dmat_scale(dmat* y,const c fac){ - dbgassert(y,NULLPTRDEREF); - for(us col=0;coln_cols;col++) { - d_scale(getdmatval(y,0,col),fac,y->n_rows); - } - -} -/** - * y = fac * y - * - * @param y - * @param fac scale factor - */ -static inline void cmat_scale(cmat* y,const c fac){ - dbgassert(y,NULLPTRDEREF); - dbgassert(y,NULLPTRDEREF); - for(us col=0;coln_cols;col++) { - c_scale(getcmatval(y,0,col),fac,y->n_rows); - } -} - -/** - * x = x + fac*y - * - * @param x - * @param y - * @param fac - */ -static inline void dmat_add_dmat(dmat* x,dmat* y,d fac) { - dbgassert(x && y,NULLPTRDEREF); - dbgassert(x->n_cols == y->n_cols,SIZEINEQUAL); - dbgassert(x->n_rows == y->n_rows,SIZEINEQUAL); - for(us col=0;coln_cols;col++) { - d_add_to(getdmatval(x,0,col), - getdmatval(y,0,col),fac,x->n_rows); - } -} -/** - * x = x + fac*y - * - * @param[in,out] x - * @param[in] y - * @param[in] fac - */ -static inline void cmat_add_cmat(cmat* x,cmat* y,c fac) { - dbgassert(x && y,NULLPTRDEREF); - dbgassert(x->n_cols == y->n_cols,SIZEINEQUAL); - dbgassert(x->n_rows == y->n_rows,SIZEINEQUAL); - for(us col=0;coln_cols;col++) { - c_add_to(getcmatval(x,0,col), - getcmatval(y,0,col),fac,x->n_rows); - } -} - -/** - * Compute the element-wise (Hadamard) product of a and b, and store - * in result - * - * @param[out] result - * @param[in] a - * @param[in] b - */ -static inline void vd_elem_prod(vd* result,const vd* a,const vd* b) { - dbgassert(result && a && b,NULLPTRDEREF); - assert_equalsize(a,b); - d_elem_prod_d(getvdval(result,0), - getvdval(a,0), - getvdval(b,0),a->n_rows); -} -/** - * Compute the element-wise (Hadamard) product of a and b, and store - * in result - * - * @param[out] result - * @param[in] a - * @param[in] b - */ -static inline void vc_hadamard(vc* result,const vc* a,const vc* b) { - fsTRACE(15); - dbgassert(result && a && b,NULLPTRDEREF); - assert_equalsize(a,b); - assert_equalsize(a,result); - - c_hadamard(getvcval(result,0), - getvcval(a,0), - getvcval(b,0), - a->n_rows); - - check_overflow_vx(*result); - check_overflow_vx(*a); - check_overflow_vx(*b); - feTRACE(15); -} - -/** - * Compute the element-wise (Hadamard) product of a and b, and store - * in result - * - * @param[out] result - * @param[in] a - * @param[in] b - */ -static inline void cmat_hadamard(cmat* result, - const cmat* a, - const cmat* b) { - fsTRACE(15); - dbgassert(result && a && b,NULLPTRDEREF); - - dbgassert(result->n_rows==a->n_rows,SIZEINEQUAL); - dbgassert(result->n_cols==a->n_cols,SIZEINEQUAL); - dbgassert(b->n_rows==a->n_rows,SIZEINEQUAL); - dbgassert(b->n_cols==a->n_cols,SIZEINEQUAL); - - for(us col=0;coln_cols;col++) { - c_hadamard(getcmatval(result,0,col), - getcmatval(a,0,col), - getcmatval(b,0,col), - a->n_rows); - } - - feTRACE(15); -} - - -/** - * Compute the matrix vector product for complex-valued types: b = A*x. - * - * @param[in] A Matrix A - * @param[in] x Vector x - * @param[out] b Result of computation - */ -void cmv_dot(const cmat* A, - const vc* restrict x, - vc* restrict b); - -int lsq_solve(const cmat* A, - const vc* restrict b, - vc* restrict x); - -/** - * Compute the norm of the difference of two complex matrices - * - * @param A - * @param B - */ -d cmat_normdiff(const cmat* A,const cmat* B); - -/** - * Computes the Kronecker product of a kron b, stores result in result. - * - * @param a a - * @param b b - * @param result a kron b - */ -void kronecker_product(const cmat* a,const cmat* b,cmat* result); - -#endif // LASP_ALG_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_alloc.h b/src/lasp/c/lasp_alloc.h deleted file mode 100644 index f6a8a55..0000000 --- a/src/lasp/c/lasp_alloc.h +++ /dev/null @@ -1,33 +0,0 @@ -// lasp_alloc.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// memory allocation functions. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_ALLOC_H -#define LASP_ALLOC_H -#include -#include "lasp_tracer.h" -/** - * Reserved words for memory allocation. Can be changed to something - * else when required. For example for debugging purposes. - */ -static inline void* a_malloc(size_t nbytes) { -#if LASP_DEBUG - if(nbytes == 0) { - FATAL("Tried to allocate 0 bytes"); - } -#endif - void* ptr = malloc(nbytes); - if(!ptr) { - FATAL("Memory allocation failed. Exiting"); - } - return ptr; -} -#define a_free free -#define a_realloc realloc - -#endif // LASP_ALLOC_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_aps.c b/src/lasp/c/lasp_aps.c deleted file mode 100644 index 75b83a5..0000000 --- a/src/lasp/c/lasp_aps.c +++ /dev/null @@ -1,291 +0,0 @@ -// lasp_aps.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-5) -#include "lasp_aps.h" -#include "lasp_ps.h" -#include "lasp_alg.h" -#include "lasp_dfifo.h" - -/* Multiplication factor for the maximum size of the fifo queue. This - * factor is multiplied by nfft to obtain the maximum size of the - * fifo. */ -#define FIFO_SIZE_MULT (5) - -typedef struct AvPowerSpectra_s { - us nfft, nchannels; - us overlap; /* Number of samples to overlap */ - - us naverages; /* Counter that counts the number of - * averages taken for the computation - * of this averaged power spectra. */ - - dmat buffer; /**< Buffer storage of some of the - * previous samples. Number of rows is - * equal to nfft. */ - - - dFifo* fifo; /* Sample fifo storage */ - - cmat* ps_storage; /**< Here we store the averaged - * results for each Cross-power - * spectra computed so far. */ - - cmat ps_result; - - cmat ps_single; /**< This is the work area for a - * PowerSpectra computation on a - * single block */ - - - vd weighting; /**< This array stores the time weighting - * coefficients for a running - * spectrogram. The vector length is - * zero for a full averaged (not - * running spectrogram). */ - - us oldest_block; /**< Index of oldest block in - * Spectrogram mode */ - - - PowerSpectra* ps; /**< Pointer to underlying - * PowerSpectra calculator. */ - -} AvPowerSpectra; - -void AvPowerSpectra_free(AvPowerSpectra* aps) { - fsTRACE(15); - - PowerSpectra_free(aps->ps); - dFifo_free(aps->fifo); - dmat_free(&aps->buffer); - - us nweight = aps->weighting.n_rows; - if(nweight > 0) { - for(us blockno = 0; blockno < nweight; blockno++) { - cmat_free(&aps->ps_storage[blockno]); - } - a_free(aps->ps_storage); - vd_free(&aps->weighting); - } - - cmat_free(&aps->ps_single); - cmat_free(&aps->ps_result); - a_free(aps); - - feTRACE(15); -} - -AvPowerSpectra* AvPowerSpectra_alloc(const us nfft, - const us nchannels, - const d overlap_percentage, - const WindowType wt, - const vd* weighting) { - - fsTRACE(15); - - /* Check nfft */ - if(nfft==0 || nfft % 2 != 0 || nfft > LASP_MAX_NFFT) { - WARN("Invalid nfft"); - feTRACE(15); - return NULL; - } - - /* Check overlap percentage */ - if(overlap_percentage < 0. || overlap_percentage >= 100.) { - WARN("Invalid overlap percentage"); - feTRACE(15); - return NULL; - } - - /* Compute and check overlap offset */ - us overlap = (us) (overlap_percentage*((d) nfft)/100); - if(overlap == nfft) { - WARN("Overlap percentage results in full overlap, decreasing overlap."); - overlap--; - } - - PowerSpectra* ps = PowerSpectra_alloc(nfft,wt); - if(!ps) { - WARN(ALLOCFAILED "ps"); - feTRACE(15); - return NULL; - } - - AvPowerSpectra* aps = a_malloc(sizeof(AvPowerSpectra)); - - aps->nchannels = nchannels; - aps->nfft = nfft; - aps->ps = ps; - aps->naverages = 0; - aps->overlap = overlap; - aps->buffer = dmat_alloc(nfft,nchannels); - aps->oldest_block = 0; - if(weighting) { - - us nweight = weighting->n_rows; - iVARTRACE(15,nweight); - /* Allocate vectors and matrices */ - aps->ps_storage = a_malloc(nweight*sizeof(cmat)); - for(us blockno = 0; blockno < nweight; blockno++) { - aps->ps_storage[blockno] = cmat_alloc(nfft/2+1,nchannels*nchannels); - cmat_set(&aps->ps_storage[blockno],0); - } - - /* Allocate space and copy weighting coefficients */ - aps->weighting = vd_alloc(weighting->n_rows); - vd_copy(&aps->weighting,weighting); - } - else { - TRACE(15,"no weighting"); - aps->weighting.n_rows = 0; - } - - aps->ps_result = cmat_alloc(nfft/2+1,nchannels*nchannels); - aps->ps_single = cmat_alloc(nfft/2+1,nchannels*nchannels); - cmat_set(&aps->ps_result,0); - - aps->fifo = dFifo_create(nchannels,FIFO_SIZE_MULT*nfft); - - feTRACE(15); - return aps; -} - - -us AvPowerSpectra_getAverages(const AvPowerSpectra* ps) { - return ps->naverages; -} - -/** - * Helper function that adds a block of time data to the APS - * - * @param aps AvPowerSpectra handle - * @param block Time data block. Size should be exactly nfft*nchannels. - */ -static void AvPowerSpectra_addBlock(AvPowerSpectra* aps, - const dmat* block) { - fsTRACE(15); - - dbgassert(aps && block,NULLPTRDEREF); - dbgassert(block->n_rows == aps->nfft,"Invalid block n_rows"); - dbgassert(block->n_cols == aps->nchannels,"Invalid block n_cols"); - - const us nfft = aps->nfft; - iVARTRACE(15,nfft); - - cmat* ps_single = &aps->ps_single; - cmat* ps_storage = aps->ps_storage; - cmat* ps_result = &aps->ps_result; - - PowerSpectra_compute(aps->ps, - block, - ps_single); - - vd weighting = aps->weighting; - us nweight = weighting.n_rows; - - if(nweight == 0) { - - /* Overall mode */ - c naverages = (++aps->naverages); - - /* Scale previous result */ - cmat_scale(ps_result, - (naverages-1)/naverages); - - /* Add new result, scaled properly */ - cmat_add_cmat(ps_result, - ps_single,1/naverages); - - } - else { - cmat_set(ps_result,0); - - us* oldest_block = &aps->oldest_block; - /* uVARTRACE(20,*oldest_block); */ - - cmat_copy(&ps_storage[*oldest_block],ps_single); - - uVARTRACE(16,*oldest_block); - if(aps->naverages < nweight) { - ++(aps->naverages); - } - - /* Update pointer to oldest block */ - (*oldest_block)++; - *oldest_block %= nweight; - - us block_index_weight = *oldest_block; - - for(us block = 0; block < aps->naverages; block++) { - /* Add new result, scaled properly */ - - c weight_fac = *getvdval(&weighting, - aps->naverages-1-block); - - /* cVARTRACE(20,weight_fac); */ - cmat_add_cmat(ps_result, - &ps_storage[block_index_weight], - weight_fac); - - block_index_weight++; - block_index_weight %= nweight; - - - } - - - - } - - feTRACE(15); -} - - -cmat* AvPowerSpectra_addTimeData(AvPowerSpectra* aps, - const dmat* timedata) { - - fsTRACE(15); - - dbgassert(aps && timedata,NULLPTRDEREF); - const us nchannels = aps->nchannels; - const us nfft = aps->nfft; - - dbgassert(timedata->n_cols == nchannels,"Invalid time data"); - dbgassert(timedata->n_rows > 0,"Invalid time data. " - "Should at least have one row"); - - /* dFifo handle */ - dFifo* fifo = aps->fifo; - dFifo_push(fifo,timedata); - - /* Temporary storage buffer */ - dmat* buffer = &aps->buffer; - - /* Pop samples from the fifo while there are still at - * least nfft samples available */ - while (dFifo_size(fifo) >= nfft) { - int popped = dFifo_pop(fifo, - buffer, - aps->overlap); /* Keep 'overlap' - * number of samples - * in the queue */ - - dbgassert((us) popped == nfft,"Bug in dFifo"); - /* Process the block of time data */ - AvPowerSpectra_addBlock(aps,buffer); - - } - - feTRACE(15); - return &aps->ps_result; -} - - - - -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_aps.h b/src/lasp/c/lasp_aps.h deleted file mode 100644 index 973037a..0000000 --- a/src/lasp/c/lasp_aps.h +++ /dev/null @@ -1,85 +0,0 @@ -// aps.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_APS_H -#define LASP_APS_H -#include "lasp_types.h" -#include "lasp_mat.h" -#include "lasp_window.h" - -typedef enum { - Linear=0, - Exponential=1 -} TimeWeighting; - - -typedef struct AvPowerSpectra_s AvPowerSpectra; - -/** - * Allocates a AvPowerSpectra object which is used to compute an - * averaged power spectra from given input time samples. - * - * @param[in] nfft Size of the fft, - * - * @param[in] nchannels Number of channels of time data to allocate - * for. - * - * @param[in] overlap_percentage. Should be 0<= overlap_pct < 100. The - * overlap percentage will be coerced, as the offset will be an - * integer number of samples at all cases. Therefore, the real overlap - * percentage can be obtained later on - * - * @param wt - * - * @return - */ -AvPowerSpectra* AvPowerSpectra_alloc(const us nfft, - const us nchannels, - const d overlap_percentage, - const WindowType wt, - const vd* spectrogram_weighting); - - -/** - * Computes the real overlap percentage, from the integer overlap - * offset. - * - * @param aps - */ -d AvPowerSpectra_realoverlap_pct(const AvPowerSpectra* aps); - -/** - * Return the current number of averages taken to obtain the result. - * - * @return The number of averages taken. - */ -us AvPowerSpectra_getAverages(const AvPowerSpectra*); - -/** - * Add time data to this Averaged Power Spectra. Remembers the part - * that should be overlapped with the next time data. - * - * @param aps AvPowerSpectra handle - * - * @param timedata Pointer to timedata buffer. Number of rows SHOULD - * *at least* be nfft. Number of columns should exactly be nchannels. - * - * @return pointer to the averaged power spectra buffer that is being - * written to. Pointer is valid until AvPowerSpectra_free() is called. - */ -cmat* AvPowerSpectra_addTimeData(AvPowerSpectra* aps, - const dmat* timedata); - -/** - * Free storage of the AvPowerSpectra - */ -void AvPowerSpectra_free(AvPowerSpectra*); - - -#endif // LASP_APS_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_assert.c b/src/lasp/c/lasp_assert.c deleted file mode 100644 index 3105e36..0000000 --- a/src/lasp/c/lasp_assert.c +++ /dev/null @@ -1,29 +0,0 @@ -// lasp_assert.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#include "lasp_config.h" -#include "lasp_assert.h" - -#ifdef LASP_DEBUG -#include -#include -#include "lasp_types.h" -#define MAX_MSG 200 - -void DBG_AssertFailedExtImplementation(const char* filename, - us linenr, - const char* extendedinfo) -{ - char scratchpad[MAX_MSG]; - - sprintf(scratchpad,"ASSERT: file %s line %lu: (%s)\n", - filename, linenr, extendedinfo); - printf("%s\n", scratchpad); - abort(); -} - -#endif diff --git a/src/lasp/c/lasp_assert.h b/src/lasp/c/lasp_assert.h deleted file mode 100644 index 2524c43..0000000 --- a/src/lasp/c/lasp_assert.h +++ /dev/null @@ -1,44 +0,0 @@ -// ascee_assert.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// Basic tools for debugging using assert statements including text. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_ASSERT_H -#define LASP_ASSERT_H - -#define OUTOFBOUNDSMATR "Out of bounds access on matrix row" -#define OUTOFBOUNDSMATC "Out of bounds access on matrix column" -#define OUTOFBOUNDSVEC "Out of bounds access on vector" -#define SIZEINEQUAL "Array sizes not equal" -#define ALLOCFAILED "Memory allocation failure in: " -#define NULLPTRDEREF "Null pointer dereference in: " - -#ifdef LASP_DEBUG -#include "lasp_types.h" - - -void DBG_AssertFailedExtImplementation(const char* file, - const us line, - const char* string); - -#define dbgassert(assertion, assert_string) \ - if (!(assertion)) \ - { \ - DBG_AssertFailedExtImplementation(__FILE__, __LINE__, assert_string ); \ - } - -#define assertvalidptr(ptr) dbgassert(ptr,NULLPTRDEREF) - -#else // LASP_DEBUG not defined - -#define dbgassert(assertion, assert_string) - -#define assertvalidptr(ptr) dbgassert(ptr,NULLPTRDEREF) - -#endif // LASP_DEBUG - -#endif // LASP_ASSERT_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_bem.c b/src/lasp/c/lasp_bem.c deleted file mode 100644 index 7d14399..0000000 --- a/src/lasp/c/lasp_bem.c +++ /dev/null @@ -1,20 +0,0 @@ -// bem.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#include "types.h" - -typedef struct { - umat elem; - vd nodes; - ElemType elemtype; -} BemMesh; - -int Ms_osc_free(BemMesh* mesh,cmat* Ms) { - - -} -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_decimation.c b/src/lasp/c/lasp_decimation.c deleted file mode 100644 index e6ab6aa..0000000 --- a/src/lasp/c/lasp_decimation.c +++ /dev/null @@ -1,221 +0,0 @@ -// lasp_decimation.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// Implementation of the decimator -////////////////////////////////////////////////////////////////////// -#include "lasp_decimation.h" -#include "lasp_firfilterbank.h" -#include "lasp_tracer.h" -#include "lasp_alloc.h" -#include "lasp_dfifo.h" - -// The Maximum number of taps in a decimation filter -#define DEC_FILTER_MAX_TAPS (128) - -// The FFT length -#define DEC_FFT_LEN (1024) - -typedef struct { - DEC_FAC df; - us dec_fac; - us ntaps; - d h[DEC_FILTER_MAX_TAPS]; -} DecFilter; - -static __thread DecFilter DecFilters[] = { - {DEC_FAC_4,4,128, -#include "dec_filter_4.c" - } -}; - - -typedef struct Decimator_s { - us nchannels; - us dec_fac; - Firfilterbank** fbs; - dFifo* output_fifo; -} Decimator; - - -Decimator* Decimator_create(us nchannels,DEC_FAC df) { - fsTRACE(15); - - /* Find the filter */ - const us nfilters = sizeof(DecFilters)/sizeof(DecFilter); - DecFilter* filter = DecFilters; - bool found = false; - for(us filterno = 0;filterno < nfilters; filterno++) { - if(filter->df == df) { - TRACE(15,"Found filter"); - found = true; - break; - } - filter++; - } - if(!found) { - WARN("Decimation factor not found in list of filters"); - return NULL; - } - - /* Create the filterbanks */ - Decimator* dec = a_malloc(sizeof(Decimator)); - 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;channelnofbs[channelno] = Firfilterbank_create(&h,DEC_FFT_LEN); - } - - dmat_free(&h); - - /* Create input and output fifo's */ - dec->output_fifo = dFifo_create(nchannels,DEC_FFT_LEN); - - feTRACE(15); - return dec; -} - -static void lasp_downsample(dmat* downsampled, - const dmat* filtered, - us dec_fac) { - fsTRACE(15); - dbgassert(downsampled && filtered,NULLPTRDEREF); - dbgassert(filtered->n_rows/dec_fac == downsampled->n_rows, - "Incompatible number of rows"); - dbgassert(downsampled->n_cols == filtered->n_cols, - "Incompatible number of rows"); - - dbgassert(dec_fac> 0,"Invalid dec_fac"); - - for(us col=0;coln_cols;col++) { - /* Low-level BLAS copy. */ - d_copy(getdmatval(downsampled,0,col), - getdmatval(filtered,0,col), - downsampled->n_rows, /* number */ - 1, /* incx out */ - dec_fac); /* incx in */ - } - - check_overflow_xmat(*downsampled); - check_overflow_xmat(*filtered); - - feTRACE(15); -} - -dmat Decimator_decimate(Decimator* dec,const dmat* samples) { - - fsTRACE(15); - dbgassert(dec && samples,NULLPTRDEREF); - const us nchannels = dec->nchannels; - const us dec_fac = dec->dec_fac; - dbgassert(samples->n_cols == nchannels,"Invalid number of " - "channels in samples"); - dbgassert(samples->n_rows > 0,"Number of rows should be >0") - - /* Not downsampled, but filtered result */ - dmat filtered; - - /* Filter each channel and store result in filtered. In first - * iteration the right size for filtered is allocated. */ - - for(us chan=0;chanfbs[chan], - &samples_channel); - - dbgassert(filtered_res.n_cols == 1,"Bug in Firfilterbank"); - - vd_free(&samples_channel); - - if(filtered_res.n_rows > 0) { - if(chan==0) { - /* Allocate space for result */ - filtered = dmat_alloc(filtered_res.n_rows, - nchannels); - - } - dmat filtered_col = dmat_submat(&filtered, - 0,chan, - filtered_res.n_rows, - 1); - - dbgassert(filtered_res.n_rows == filtered_col.n_rows, - "Not all Firfilterbank's have same output number" - " of rows!"); - - dmat_copy_rows(&filtered_col, - &filtered_res, - 0,0,filtered_res.n_rows); - - dmat_free(&filtered_res); - dmat_free(&filtered_col); - } - else { - filtered = dmat_alloc(0, nchannels); - } - } - if(filtered.n_rows > 0) { - /* Push filtered result into output fifo */ - dFifo_push(dec->output_fifo, - &filtered); - } - dmat_free(&filtered); - - - /* Now, downsample stuff */ - dmat downsampled; - uVARTRACE(15,dec_fac); - us fifo_size = dFifo_size(dec->output_fifo); - if((fifo_size / dec_fac) > 0) { - - filtered = dmat_alloc((fifo_size/dec_fac)*dec_fac, - nchannels); - - int nsamples = dFifo_pop(dec->output_fifo, - &filtered,0); - - dbgassert((us) nsamples % dec_fac == 0 && nsamples > 0, - "BUG"); - dbgassert((us) nsamples == (fifo_size/dec_fac)*dec_fac,"BUG"); - - downsampled = dmat_alloc(nsamples/dec_fac,nchannels); - /* Do the downsampling work */ - lasp_downsample(&downsampled,&filtered,dec_fac); - - dmat_free(&filtered); - } - else { - TRACE(15,"Empty downsampled"); - downsampled = dmat_alloc(0,0); - } - - feTRACE(15); - /* return filtered; */ - return downsampled; -} - - -void Decimator_free(Decimator* dec) { - fsTRACE(15); - dbgassert(dec,NULLPTRDEREF); - dFifo_free(dec->output_fifo); - - for(us chan=0;channchannels;chan++) { - Firfilterbank_free(dec->fbs[chan]); - } - - a_free(dec->fbs); - a_free(dec); - - feTRACE(15); -} -////////////////////////////////////////////////////////////////////// - diff --git a/src/lasp/c/lasp_decimation.h b/src/lasp/c/lasp_decimation.h deleted file mode 100644 index 1115502..0000000 --- a/src/lasp/c/lasp_decimation.h +++ /dev/null @@ -1,54 +0,0 @@ -// lasp_decimation.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: Sample decimator, works on a whole block of -// (uninterleaved) time samples at once. Decimates (downsamples and -// low-pass filters by a factor given integer factor. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_DECIMATION_H -#define LASP_DECIMATION_H -#include "lasp_types.h" -#include "lasp_mat.h" - -typedef struct Decimator_s Decimator; - -typedef enum DEC_FAC_e { - DEC_FAC_4 = 0, // Decimate by a factor of 4 -} DEC_FAC; - -/** - * Create a decimation filter for a given number of channels and - * decimation factor - * - * @param nchannels Number of channels - - * @param d Decimation factor. Should be one of the implemented - * ones. (TODO: add more. At this point we have only a decimation - * factor of 4 implemented) - * - * @return Decimator handle. NULL on error - */ -Decimator* Decimator_create(us nchannels,DEC_FAC d); - -/** - * Decimate given samples - * - * @param samples - * - * @return Decimated samples. Can be an empty array. - */ -dmat Decimator_decimate(Decimator* dec,const dmat* samples); - - -d Decimator_get_cutoff(Decimator*); -/** - * Free memory corresponding to Decimator - * - * @param dec Decimator handle. - */ -void Decimator_free(Decimator* dec); - -#endif // LASP_DECIMATION_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_dfifo.c b/src/lasp/c/lasp_dfifo.c deleted file mode 100644 index d5e2e0e..0000000 --- a/src/lasp/c/lasp_dfifo.c +++ /dev/null @@ -1,145 +0,0 @@ -// lasp_dfifo.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// Implementation of the dFifo queue -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-10) -#include "lasp_dfifo.h" - -typedef struct dFifo_s { - dmat queue; - us start_row; - us end_row; -} dFifo; -us dFifo_size(dFifo* fifo) { - fsTRACE(15); - dbgassert(fifo,NULLPTRDEREF); - dbgassert(fifo->start_row <= fifo->end_row,"BUG"); - feTRACE(15); - return fifo->end_row-fifo->start_row; -} - -/** - * Change the max size of the dFifo to a new max size specified. Max size - * should be larger than fifo size. Resets start row to 0 - * - * @param fifo - * @param new_size - */ -static void dFifo_change_maxsize(dFifo* fifo,const us new_max_size) { - fsTRACE(15); - dmat old_queue = fifo->queue; - - dbgassert(new_max_size >= dFifo_size(fifo),"BUG"); - const us size = dFifo_size(fifo); - - dmat new_queue = dmat_alloc(new_max_size,old_queue.n_cols); - if(size > 0) { - dmat_copy_rows(&new_queue, - &old_queue, - 0, - fifo->start_row, - size); - } - - dmat_free(&old_queue); - fifo->queue = new_queue; - fifo->end_row -= fifo->start_row; - fifo->start_row = 0; - - feTRACE(15); -} - -dFifo* dFifo_create(const us nchannels, - const us init_size) { - - fsTRACE(15); - dFifo* fifo = a_malloc(sizeof(dFifo)); - fifo->queue = dmat_alloc(init_size,nchannels); - fifo->start_row = 0; - fifo->end_row = 0; - feTRACE(15); - return fifo; -} -void dFifo_free(dFifo* fifo) { - fsTRACE(15); - dmat_free(&fifo->queue); - a_free(fifo); - feTRACE(15); -} -void dFifo_push(dFifo* fifo,const dmat* data) { - fsTRACE(15); - dbgassert(fifo && data, NULLPTRDEREF); - dbgassert(data->n_cols == fifo->queue.n_cols, - "Invalid number of columns in data"); - - - const us added_size = data->n_rows; - - dmat queue = fifo->queue; - const us max_size = queue.n_rows; - - us* end_row = &fifo->end_row; - - if(added_size + dFifo_size(fifo) > max_size) { - dFifo_change_maxsize(fifo,2*(max_size+added_size)); - - /* Now the stack of this function is not valid anymore. Best - * thing to do is restart the function. */ - dFifo_push(fifo,data); - feTRACE(15); - return; - } - else if(*end_row + added_size > max_size) { - dFifo_change_maxsize(fifo,max_size); - /* Now the stack of this function is not valid anymore. Best - * thing to do is restart the function. */ - dFifo_push(fifo,data); - feTRACE(15); - return; - - } - - /* Now, copy samples */ - dmat_copy_rows(&queue, /* to */ - data, /* from */ - *end_row, /* startrow_to */ - 0, /* startrow_from */ - added_size); /* n_rows */ - - /* Increase the size */ - *end_row += added_size; - - feTRACE(15); -} -int dFifo_pop(dFifo* fifo,dmat* data,const us keep) { - fsTRACE(15); - dbgassert(fifo && data,NULLPTRDEREF); - dbgassert(data->n_cols == fifo->queue.n_cols, - "Invalid number of columns in data"); - dbgassert(keep < data->n_rows, "Number of samples to keep should" - " be smaller than requested number of samples"); - - us* start_row = &fifo->start_row; - us cur_size = dFifo_size(fifo); - us requested = data->n_rows; - - us obtained = requested > cur_size ? cur_size : requested; - dbgassert(obtained > keep,"Number of samples to keep should be" - " smaller than requested number of samples"); - - - dmat_copy_rows(data, - &fifo->queue, - 0, - *start_row, - obtained); - - *start_row += obtained - keep; - - feTRACE(15); - return (int) obtained; -} -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_dfifo.h b/src/lasp/c/lasp_dfifo.h deleted file mode 100644 index 98f198f..0000000 --- a/src/lasp/c/lasp_dfifo.h +++ /dev/null @@ -1,73 +0,0 @@ -// lasp_dfifo.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// API of a contiguous fifo buffer of samples. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_DFIFO_H -#define LASPDFIFO_H -#include "lasp_types.h" -#include "lasp_mat.h" - -typedef struct dFifo_s dFifo; - -/** - * Create a fifo buffer - * - * @param nchannels Number of channels to store for - * @param max_size Maximum size of the queue. - * - * @return Pointer to fifo queue (no null ptr) - */ -dFifo* dFifo_create(const us nchannels, - const us init_size); - - -/** - * Pushes samples into the fifo. - * - * @param fifo dFifo handle - * - * @param data data to push. Number of columns should be equal to - * nchannels. - */ -void dFifo_push(dFifo* fifo,const dmat* data); - -/** - * Pop samples from the queue - * - * @param[in] fifo dFifo handle - - * @param[out] data Pointer to dmat where popped data will be - * stored. Should have nchannels number of columns. If n_rows is - * larger than current storage, the queue is emptied. - - * @param[in] keep Keeps a number of samples for the next dFifo_pop(). If - * keep=0, then no samples will be left. Keep should be smaller than - * the number of rows in data. - * - * @return Number of samples obtained in data. - */ -int dFifo_pop(dFifo* fifo,dmat* data,const us keep); - -/** - * Returns current size of the fifo - * - * @param[in] fifo dFifo handle - * - * @return Current size - */ -us dFifo_size(dFifo* fifo); - -/** - * Free a dFifo object - * - * @param[in] fifo dFifo handle. - */ -void dFifo_free(dFifo* fifo); - -#endif // LASP_DFIFO_H -////////////////////////////////////////////////////////////////////// - diff --git a/src/lasp/c/lasp_eq.c b/src/lasp/c/lasp_eq.c deleted file mode 100644 index 1cf462f..0000000 --- a/src/lasp/c/lasp_eq.c +++ /dev/null @@ -1,70 +0,0 @@ -#include "lasp_eq.h" -#include "lasp_assert.h" - -typedef struct Eq { - Sosfilterbank* fb; - us nfilters; - vd ampl_values; -} Eq; - -Eq* Eq_create(Sosfilterbank* fb) { - fsTRACE(15); - assertvalidptr(fb); - Eq* eq = a_malloc(sizeof(Eq)); - eq->fb = fb; - eq->nfilters = Sosfilterbank_getFilterbankSize(fb); - eq->ampl_values = vd_alloc(eq->nfilters); - vd_set(&(eq->ampl_values), 1.0); - feTRACE(15); - return eq; -} -vd Eq_equalize(Eq* eq,const vd* input_data) { - fsTRACE(15); - assertvalidptr(eq); - assert_vx(input_data); - vd result = vd_alloc(input_data->n_rows); - dmat_set(&result, 0); - dmat filtered = Sosfilterbank_filter(eq->fb, input_data); - - for(us filter=0;filternfilters;filter++) { - d ampl = *getvdval(&(eq->ampl_values), filter); - /// TODO: Replace this code with something more fast from BLAS. - for(us sample=0;samplen_rows == eq->nfilters, "Invalid levels size"); - for(us ch=0;chnfilters;ch++){ - d level = *getvdval(levels, ch); - *getvdval(&(eq->ampl_values), ch) = d_pow(10, level/20); - } - - feTRACE(15); -} - -us Eq_getNLevels(const Eq* eq) { - fsTRACE(15); - assertvalidptr(eq); - feTRACE(15); - return eq->nfilters; -} - -void Eq_free(Eq* eq) { - fsTRACE(15); - assertvalidptr(eq); - assertvalidptr(eq->fb); - Sosfilterbank_free(eq->fb); - vd_free(&(eq->ampl_values)); - feTRACE(15); -} diff --git a/src/lasp/c/lasp_eq.h b/src/lasp/c/lasp_eq.h deleted file mode 100644 index 1afefff..0000000 --- a/src/lasp/c/lasp_eq.h +++ /dev/null @@ -1,59 +0,0 @@ -// lasp_eq.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: Implementation of an equalizer using the Second Order Sections -// filter bank implementation. Applies all filterbanks in parallel and -// recombines to create a single output signal. -#pragma once -#ifndef LASP_EQ_H -#define LASP_EQ_H -#include "lasp_sosfilterbank.h" -typedef struct Eq Eq; - -/** - * Initialize an equalizer using the given Filterbank. Note: takes over pointer - * ownership of fb! Sets levels of all filterbanks initially to 0 dB. - * - * @param[in] fb Filterbank to be used in equalizer - * @return Equalizer handle, NULL on error - * */ -Eq* Eq_create(Sosfilterbank* fb); - -/** - * Equalize a given piece of data using current settings of the levels. - * - * @param[in] eq Equalizer handle - * @param[in] input_data Input data to equalize - * @return Equalized data. Newly allocated vector. Ownership is transferred. - * */ -vd Eq_equalize(Eq* eq,const vd* input_data); - -/** - * Returns number of channels of the equalizer. Note: takes over pointer - * ownership of fb! - * - * @param[in] eq Equalizer handle - * */ -us Eq_getNLevels(const Eq* eq); - -/** - * Set amplification values for each filter in the equalizer. - * - * @param[in] eq Equalizer handle - * @param[in] levels: Vector with level values for each channel. Should have - * length equal to the number of filters in the filterbank. - * */ -void Eq_setLevels(Eq* eq, const vd* levels); - -/** - * Cleans up an existing Equalizer - * - * @param[in] eq Equalizer handle - */ -void Eq_free(Eq* eq); - - -#endif // LASP_EQ_H -// ////////////////////////////////////////////////////////////////////// - diff --git a/src/lasp/c/lasp_fft.h b/src/lasp/c/lasp_fft.h deleted file mode 100644 index 5b1f39d..0000000 --- a/src/lasp/c/lasp_fft.h +++ /dev/null @@ -1,109 +0,0 @@ -// lasp_fft.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// Interface to the FFT library, multiple channel FFT's -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_FFT_H -#define LASP_FFT_H -#include "lasp_types.h" -#include "lasp_mat.h" -#ifdef __cplusplus -extern "C" { -#endif - -/** - * Load wisdom data from a NULL-terminated character string. Note that at this - * point in time this is only used by the FFTW fft backend. - * - * @param[in] wisdom: Wisdom string - */ -void load_fft_wisdom(const char* wisdom); - -/** - * Creates a NULL-terminated string containing FFTW wisdom, or state knowledge - * from other libraries. - * - * @returns A NULL-terminated string. *Should be deallocated after being used* - */ -char* store_fft_wisdom(void); - - -/** - * Perform forward FFT's on real time data. - * - */ -typedef struct Fft_s Fft; - -/** - * Construct an Fft object - * - * @param nfft Nfft size - * - * @return Pointer to Fft handle, NULL on error - */ -Fft* Fft_create(const us nfft); -/** - * Returns the nfft for this Fft instance - * - * @return nfft - */ -us Fft_nfft(const Fft*); - -/** - * Compute the fft of a single channel. - * - * @param[in] fft Fft handle. - * - * @param[in] timedata Input time data pointer, should have size nfft - * @param[out] result Pointer to result vector should have size - * nfft/2+1 - */ -void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result); - -/** - * Compute the fft of the data matrix, first axis is assumed to be - * the time axis. - * - * @param[in] fft Fft handle. - - * @param[in] timedata Input time data. First axis is assumed to be - * the time, second the channel number. - * - * @param[out] result: Fft't data, should have size (nfft/2+1) * - * nchannels - */ -void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result); - -/** - * Perform inverse fft on a single channel. - * - * @param[in] fft Fft handle. - * @param[in] freqdata Frequency domain input data, to be iFft'th. - * @param[out] timedata: iFft't data, should have size (nfft). - */ -void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* timedata); - -/** - * Perform inverse FFT - * - * @param[in] fft Fft handle - * @param[in] freqdata Frequency domain data - * @param[out] timedata Time domain result - */ -void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata); - -/** - * Free up resources of Fft handle. - * - * @param fft Fft handle. - */ -void Fft_free(Fft* fft); - -#ifdef __cplusplus -} -#endif -#endif // LASP_FFT_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_firfilterbank.c b/src/lasp/c/lasp_firfilterbank.c deleted file mode 100644 index ad3bebe..0000000 --- a/src/lasp/c/lasp_firfilterbank.c +++ /dev/null @@ -1,174 +0,0 @@ -// lasp_filterbank.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// Firfilterbank implementation. -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-5) -#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 Firfilterbank_s { - us nfft; - - us P_m_1; /**< Filter length minus one */ - - cmat filters; /* Frequency-domain filter - * coefficients */ - dFifo* output_fifo; - dFifo* input_fifo; - - Fft* fft; /* Handle to internal FFT-function */ - -} Firfilterbank; - -Firfilterbank* Firfilterbank_create(const dmat* h, - const us nfft) { - - fsTRACE(15); - dbgassert(h,NULLPTRDEREF); - const us P = h->n_rows; - const us nfilters = h->n_cols; - - if(P > nfft/2) { - WARN("Filter order should be <= nfft/2"); - return NULL; - } - - Fft* fft = Fft_create(nfft); - if(!fft) { - WARN("Fft allocation failed"); - return NULL; - } - - Firfilterbank* fb = a_malloc(sizeof(Firfilterbank)); - - fb->nfft = nfft; - fb->P_m_1 = P-1; - fb->fft = fft; - fb->filters = cmat_alloc(nfft/2+1,nfilters); - - fb->output_fifo = dFifo_create(nfilters,FIFO_SIZE_MULT*nfft); - fb->input_fifo = dFifo_create(1,FIFO_SIZE_MULT*nfft); - - // Initialize the input fifo with zeros. - // dmat init_zero = dmat_alloc(nfft - P, 1); - // dmat_set(&init_zero,0); - // dFifo_push(fb->input_fifo, &init_zero); - // dmat_free(&init_zero); - - /* Create a temporary buffer which is going to be FFT'th to - * contain the filter transfer functions. - */ - dmat temp = dmat_alloc(nfft,nfilters); - dmat_set(&temp,0); - dmat_copy_rows(&temp,h,0,0,h->n_rows); - - /* Fft the FIR impulse responses */ - Fft_fft(fb->fft,&temp,&fb->filters); - - dmat_free(&temp); - - feTRACE(15); - return fb; -} -void Firfilterbank_free(Firfilterbank* fb) { - fsTRACE(15); - dbgassert(fb,NULLPTRDEREF); - cmat_free(&fb->filters); - dFifo_free(fb->input_fifo); - dFifo_free(fb->output_fifo); - Fft_free(fb->fft); - a_free(fb); - feTRACE(15); -} -dmat Firfilterbank_filter(Firfilterbank* fb, - const vd* x) { - - fsTRACE(15); - dbgassert(fb && x ,NULLPTRDEREF); - - dFifo* input_fifo = fb->input_fifo; - dFifo* output_fifo = fb->output_fifo; - - const us nfft = fb->nfft; - const us nfilters = fb->filters.n_cols; - - /* Push samples to the input fifo */ - dFifo_push(fb->input_fifo,x); - - dmat input_block = dmat_alloc(nfft,1); - - /* FFT'th filter coefficients */ - cmat input_fft = cmat_alloc(nfft/2+1,1); - - /* Output of the fast convolution */ - cmat output_fft = cmat_alloc(nfft/2+1,nfilters); - - /* Inverse FFT'th output */ - dmat output_block = dmat_alloc(nfft,nfilters); - - while (dFifo_size(input_fifo) >= nfft) { - - us nsamples = dFifo_pop(input_fifo, - &input_block, - fb->P_m_1 /* save P-1 samples */ - ); - - dbgassert(nsamples == nfft,"BUG in dFifo"); - - Fft_fft(fb->fft,&input_block,&input_fft); - - vc input_fft_col = cmat_column(&input_fft,0); - - for(us col=0;colfilters,col); - - vc_hadamard(&output_fft_col, - &input_fft_col, - &filter_col); - - vc_free(&output_fft_col); - vc_free(&filter_col); - - } - - vc_free(&input_fft_col); - - Fft_ifft(fb->fft,&output_fft,&output_block); - - dmat valid_samples = dmat_submat(&output_block, - fb->P_m_1,0, /* startrow, startcol */ - nfft-fb->P_m_1, /* Number of rows */ - output_block.n_cols); - - /* Push the valid samples to the output FIFO */ - dFifo_push(fb->output_fifo,&valid_samples); - dmat_free(&valid_samples); - - } - - dmat_free(&input_block); - cmat_free(&input_fft); - cmat_free(&output_fft); - dmat_free(&output_block); - - us samples_done = dFifo_size(output_fifo); - uVARTRACE(15,samples_done); - dmat filtered_result = dmat_alloc(samples_done,nfilters); - if(samples_done) { - us samples_done2 = dFifo_pop(output_fifo,&filtered_result,0); - dbgassert(samples_done2 == samples_done,"BUG in dFifo"); - } - feTRACE(15); - return filtered_result; -} - - -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_firfilterbank.h b/src/lasp/c/lasp_firfilterbank.h deleted file mode 100644 index e0f646b..0000000 --- a/src/lasp/c/lasp_firfilterbank.h +++ /dev/null @@ -1,59 +0,0 @@ -// lasp_firfilterbank.h -// -// Author: J.A. de Jong - ASCEE -// -// 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 -// performed only over the input data, and the backwards transfer for -// each filter in the filterbank. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_FIRFILTERBANK_H -#define LASP_FIRFILTERBANK_H -#include "lasp_types.h" -#include "lasp_mat.h" -typedef struct Firfilterbank_s Firfilterbank; - -/** - * 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. - * - * @param h: matrix with filter coefficients of each of the - * filters. First axis is the axis of the filter coefficients, second - * axis is the filter number. Maximum length of the filter is nfft/2. - * - * @param nfft: FTT length for fast convolution. For good performance, - * nfft should be chosen as the nearest power of 2, approximately four - * times the filter lengths. For the lowest possible latency, it is - * better to set nfft at twice the filter length. - * - * @return Firfilterbank handle, NULL on error. - */ -Firfilterbank* Firfilterbank_create(const dmat* h,const us nfft); - -/** - * 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 Firfilterbank_filter(Firfilterbank* fb, - const vd* x); - -/** - * Cleans up an existing filter bank. - * - * @param f Filter handle - */ -void Firfilterbank_free(Firfilterbank* f); - - -#endif // LASP_FIRFILTERBANK_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_mat.c b/src/lasp/c/lasp_mat.c deleted file mode 100644 index bece9d3..0000000 --- a/src/lasp/c/lasp_mat.c +++ /dev/null @@ -1,54 +0,0 @@ -// lasp_mat.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-10) -#include "lasp_mat.h" -#include "lasp_assert.h" -#include "lasp_tracer.h" - -#include - -#if LASP_DEBUG == 1 -void print_dmat(const dmat* m) { - fsTRACE(50); - size_t row,col; - for(row=0;rown_rows;row++){ - indent_trace(); - for(col=0;coln_cols;col++){ - d val = *getdmatval(m,row,col); - printf("%c%2.2e ", val<0?'-':' ' ,d_abs(val)); - - } - printf("\n"); - - } - feTRACE(50); -} -void print_cmat(const cmat* m) { - fsTRACE(50); - size_t row,col; - for(row=0;rown_rows;row++){ - indent_trace(); - for(col=0;coln_cols;col++){ - c val = *getcmatval(m,row,col); - - d rval = creal(val); - d ival = cimag(val); - - printf("%c%2.2e%c%2.2ei ",rval< 0 ?'-': ' ', - d_abs(rval),ival<0 ? '-' : '+',d_abs(ival) ) ; - - } - printf("\n"); - - } - feTRACE(50); -} - -#endif // LASP_DEBUG == 1 - -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_mat.h b/src/lasp/c/lasp_mat.h deleted file mode 100644 index 3a349ae..0000000 --- a/src/lasp/c/lasp_mat.h +++ /dev/null @@ -1,597 +0,0 @@ -// lasp_mat.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: Basic routines for allocating, setting, freeing and -// copying of matrices and vectors. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_MAT_H -#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" - -/// Dense matrix of floating point values -typedef struct { - us n_rows; - us n_cols; - bool _foreign_data; - us colstride; - d* _data; -} dmat; - -/// Dense matrix of complex floating point values -typedef struct { - us n_rows; - us n_cols; - bool _foreign_data; - us colstride; - c* _data; -} cmat; - -typedef dmat vd; -typedef cmat vc; - -#define assert_equalsize(a,b) \ - dbgassert((a)->n_rows == (b)->n_rows,SIZEINEQUAL); \ - dbgassert((a)->n_cols == (b)->n_cols,SIZEINEQUAL); - -#define is_vx(vx) ((vx)->n_cols == 1) -#define assert_vx(vx) dbgassert(is_vx(vx),"Not a vector!") - -#define setvecval(vec,index,val) \ - assert_vx(vec); \ - dbgassert((((us) index) < (vec)->n_rows),OUTOFBOUNDSVEC); \ - (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)->colstride+(row)] = val; - -/** - * Return pointer to a value from a vector - * - * @param mat The vector - * @param row The row - */ -/** - * Return a value from a matrix of floating points - * - * @param mat The matrix - * @param row The row - * @param col The column - */ -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->colstride+row]; -} -/** - * Return a value from a matrix of complex floating points - * - * @param mat The matrix - * @param row The row - * @param col The column - */ -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->colstride+row]; -} - -static inline d* getvdval(const vd* vec,us row){ - dbgassert(vec,NULLPTRDEREF); - assert_vx(vec); - return getdmatval(vec,row,0); -} - -/** - * Return pointer to a value from a complex vector - * - * @param mat The vector - * @param row The row - */ -static inline c* getvcval(const vc* vec,us row){ - dbgassert(vec,NULLPTRDEREF); - assert_vx(vec); - return getcmatval(vec,row,0); -} - - - -#ifdef LASP_DEBUG -#define OVERFLOW_MAGIC_NUMBER (-10e-45) - -#define check_overflow_xmat(xmat) \ - TRACE(15,"Checking overflow " #xmat); \ - if(!(xmat)._foreign_data) { \ - dbgassert((xmat)._data[((xmat).n_cols-1)*(xmat).colstride+(xmat).n_rows] \ - == OVERFLOW_MAGIC_NUMBER, \ - "Buffer overflow detected on" #xmat ); \ - } \ - -#define check_overflow_vx check_overflow_xmat - -#else -#define check_overflow_vx(vx) -#define check_overflow_xmat(xmat) -#endif - -/** - * Sets all values in a matrix to the value - * - * @param mat The matrix to set - * @param value - */ -static inline void dmat_set(dmat* mat,const d value){ - dbgassert(mat,NULLPTRDEREF); - if(islikely(mat->n_cols * mat->n_rows > 0)) { - for(us col=0;coln_cols;col++) { - d_set(getdmatval(mat,0,col),value,mat->n_rows); - } - } -} -#define vd_set dmat_set - -/** - * Sets all values in a matrix to the value - * - * @param mat The matrix to set - * @param value - */ -static inline void cmat_set(cmat* mat,const c value){ - dbgassert(mat,NULLPTRDEREF); - if(islikely(mat->n_cols * mat->n_rows > 0)) { - for(us col=0;coln_cols;col++) { - c_set(getcmatval(mat,0,col),value,mat->n_rows); - } - } -} -#define vc_set cmat_set - -/** - * Allocate data for a matrix of floating points - * - * @param n_rows Number of rows - * @param n_cols Number of columns - * @param p Memory pool - * - * @return dmat with allocated data - */ -static inline dmat dmat_alloc(us n_rows, - us n_cols) { - 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)); - result._data[n_rows*n_cols] = OVERFLOW_MAGIC_NUMBER; - #else - result._data = (d*) a_malloc((n_rows*n_cols)*sizeof(d)); - #endif // LASP_DEBUG - - #ifdef LASP_DEBUG - dmat_set(&result,NAN); - #endif // LASP_DEBUG - - return result; -} - -/** - * Allocate a matrix of complex floating points - * - * @param n_rows Number of rows - * @param n_cols Number of columns - * @param p Memory pool - * - * @return cmat with allocated data - */ -static inline cmat cmat_alloc(const us n_rows, - const us n_cols) { - cmat result = { n_rows, n_cols, false, n_rows, NULL}; - - #ifdef LASP_DEBUG - result._data = (c*) a_malloc((n_rows*n_cols+1)*sizeof(c)); - result._data[n_rows*n_cols] = OVERFLOW_MAGIC_NUMBER; - #else - result._data = (c*) a_malloc((n_rows*n_cols)*sizeof(c)); - #endif // LASP_DEBUG - - #ifdef LASP_DEBUG - cmat_set(&result,NAN+I*NAN); - #endif // LASP_DEBUG - return result; -} - -/** - * Allocate data for a float vector. - * - * @param size Size of the vector - * - * @return vd with allocated data - */ -static inline vd vd_alloc(us size) { - return dmat_alloc(size,1); -} - -/** - * Allocate data for a complex vector. - * - * @param size Size of the vector - * - * @return vc with allocated data - */ -static inline vc vc_alloc(us size) { - return cmat_alloc(size,1); -} - - -/** - * 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 colstride equals to n_rows. - * - * @param n_rows Number of rows - * @param n_cols Number of columns - * @param data - * - * @return - */ -static inline dmat dmat_foreign(dmat* other) { - dbgassert(other,NULLPTRDEREF); - dmat result = {other->n_rows, - other->n_cols, - true, - other->colstride, - other->_data}; - return result; -} -/** - * Create a dmat from foreign data. Assumes the colstride of the data is - * n_rows. - * - * @param n_rows Number of rows - * @param n_cols Number of columns - * @param data Pointer to data storage - * - * @return dmat - */ -static inline dmat dmat_foreign_data(us n_rows, - us n_cols, - d* data, - bool own_data) { - - dbgassert(data,NULLPTRDEREF); - dmat result = {n_rows, - n_cols, - !own_data, - n_rows, - data}; - return result; -} -/** - * Create a cmat from foreign data. Assumes the colstride of the data is - * n_rows. - * - * @param n_rows Number of rows - * @param n_cols Number of columns - * @param data Pointer to data storage - * - * @return dmat - */ -static inline cmat cmat_foreign_data(us n_rows, - us n_cols, - c* data, - bool own_data) { - - dbgassert(data,NULLPTRDEREF); - cmat result = {n_rows, - n_cols, - !own_data, - n_rows, - data}; - return result; -} - -/** - * 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 colstride equals to n_rows. - * - * @param n_rows - * @param n_cols - * @param data - * - * @return - */ -static inline cmat cmat_foreign(cmat* other) { - dbgassert(other,NULLPTRDEREF); - cmat result = {other->n_rows, - other->n_cols, - true, - other->colstride, - other->_data}; - return result; -} - - - -/** - * Free's data of dmat. Safe to run on sub-matrices as well. - * - * @param m Matrix to free - */ -static inline void dmat_free(dmat* m) { - dbgassert(m,NULLPTRDEREF); - if(!(m->_foreign_data)) a_free(m->_data); -} -#define vd_free dmat_free - -/** - * Free's data of dmat. Safe to run on sub-matrices as well. - * - * @param m Matrix to free - */ -static inline void cmat_free(cmat* m) { - dbgassert(m,NULLPTRDEREF); - if(!(m->_foreign_data)) a_free(m->_data); -} -#define vc_free cmat_free - -/** - * Copy some rows from one matrix to another - * - * @param to Matrix to copy to - * @param from Matrix to copy from - * @param startrow_from Starting row where to get the values - * @param startrow_to Starting row where to insert the values - * @param nrows Number of rows to copy - */ -static inline void dmat_copy_rows(dmat* to,const dmat* from, - us startrow_to, - us startrow_from, - us nrows) { - us col,ncols = to->n_cols; - dbgassert(to && from,NULLPTRDEREF); - dbgassert(to->n_cols == from->n_cols,SIZEINEQUAL); - dbgassert(startrow_from+nrows <= from->n_rows,OUTOFBOUNDSMATR); - dbgassert(startrow_to+nrows <= to->n_rows,OUTOFBOUNDSMATR); - - for(col=0;coln_rows,OUTOFBOUNDSMATR); - dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC); - - dmat result = { n_rows,n_cols, - true, // Foreign data = true - parent->n_rows, // This is the colstride to get to - // the next column. - getdmatval(parent,startrow,startcol)}; - - return result; -} -/** - * Allocate a sub-matrix view of the parent - * - * @param parent Parent matrix - * @param startrow Startrow - * @param startcol Start column - * @param n_rows Number of rows in sub-matrix - * @param n_cols Number of columns in sub-matrix - * - * @return submatrix view - */ -static inline cmat cmat_submat(cmat* parent, - const us startrow, - const us startcol, - const us n_rows, - const us n_cols) { - dbgassert(false,"untested"); - dbgassert(parent,NULLPTRDEREF); - dbgassert(n_rows+startrow <= parent->n_rows,OUTOFBOUNDSMATR); - dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC); - - - cmat result = { n_rows,n_cols, - true, // Foreign data = true - parent->n_rows, // This is the colstride to get to - // the next column. - getcmatval(parent,startrow,startcol)}; - - return result; -} - -/** - * Copy contents of one matrix to another. Sizes should be equal - * - * @param to - * @param from - */ -static inline void dmat_copy(dmat* to,const dmat* from) { - dbgassert(to && from,NULLPTRDEREF); - dbgassert(to->n_rows==from->n_rows,SIZEINEQUAL); - dbgassert(to->n_cols==from->n_cols,SIZEINEQUAL); - for(us col=0;coln_cols;col++) { - d_copy(getdmatval(to,0,col), - getdmatval(from,0,col), - to->n_rows,1,1); - } -} - -/** - * Allocate a new array, with size based on other. - * - * @param[in] from: Array to copy - */ -static inline dmat dmat_alloc_from_dmat(const dmat* from) { - assertvalidptr(from); - dmat thecopy = dmat_alloc(from->n_rows, from->n_cols); - return thecopy; -} - -/** - * Copy contents of one matrix to another. Sizes should be equal - * - * @param to - * @param from - */ -static inline void cmat_copy(cmat* to,const cmat* from) { - dbgassert(to && from,NULLPTRDEREF); - dbgassert(to->n_rows==from->n_rows,SIZEINEQUAL); - dbgassert(to->n_cols==from->n_cols,SIZEINEQUAL); - for(us col=0;coln_cols;col++) { - c_copy(getcmatval(to,0,col), - getcmatval(from,0,col), - to->n_rows); - } -} - -/** - * Copy contents of one vector to another - * - * @param to : Vector to write to - * @param from : Vector to read from - */ -static inline void vd_copy(vd* to,const vd* from) { - dbgassert(to && from,NULLPTRDEREF); - assert_vx(to); - assert_vx(from); - dmat_copy(to,from); -} -/** - * Copy contents of one vector to another - * - * @param to : Vector to write to - * @param from : Vector to read from - */ -static inline void vc_copy(vc* to,const vc* from) { - dbgassert(to && from,NULLPTRDEREF); - assert_vx(to); - assert_vx(from); - cmat_copy(to,from); -} - - -/** - * Get a reference to a column of a matrix as a vector - * - * @param x Matrix - * @param col Column number - * - * @return vector with reference to column - */ -static inline vd dmat_column(dmat* x,us col) { - vd res = {x->n_rows,1,true,x->colstride,getdmatval(x,0,col)}; - return res; -} - -/** - * Get a reference to a column of a matrix as a vector - * - * @param x Matrix - * @param col Column number - * - * @return vector with reference to column - */ -static inline vc cmat_column(cmat* x,us col) { - vc res = {x->n_rows,1,true,x->colstride,getcmatval(x,0,col)}; - return res; -} - -/** - * Compute the complex conjugate of b and store result in a - * - * @param a - * @param b - */ -static inline void cmat_conj(cmat* a,const cmat* b) { - fsTRACE(15); - dbgassert(a && b,NULLPTRDEREF); - dbgassert(a->n_cols == b->n_cols,SIZEINEQUAL); - dbgassert(a->n_rows == b->n_rows,SIZEINEQUAL); - for(us col=0;coln_cols;col++) { - carray_conj(getcmatval(a,0,col),getcmatval(b,0,col),a->n_rows); - } - feTRACE(15); -} - -/** - * Take the complex conjugate of x, in place - * - * @param x - */ -static inline void cmat_conj_inplace(cmat* x) { - dbgassert(x,NULLPTRDEREF); - for(us col=0;coln_cols;col++) { - c_conj_inplace(getcmatval(x,0,col),x->n_rows); - } -} - -/** - * Computes the maximum value for each row, returns a vector with maximum - * values for each column in the matrix. - * - * @param x - */ -static inline vd dmat_max(const dmat x) { - vd max_vals = vd_alloc(x.n_cols); - d max_val = -d_inf; - for(us j=0; j< x.n_cols; j++) { - for(us i=0; i< x.n_rows; i++) { - max_val = *getdmatval(&x, i, j) < max_val? max_val : *getdmatval(&x, i,j); - } - *getvdval(&max_vals, j) = max_val; - } - return max_vals; -} - - -#if LASP_DEBUG == 1 - -void print_cmat(const cmat* m); -void print_dmat(const dmat* m); -#define print_vc(x) assert_vx(x) print_cmat(x) -#define print_vd(x) assert_vx(x) print_dmat(x) - -#else -#define print_cmat(m) -#define print_dmat(m) -#define print_vc(m) -#define print_vd(m) - -#endif // LASP_DEBUG == 1 - -#endif // LASP_MAT_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_math_raw.c b/src/lasp/c/lasp_math_raw.c deleted file mode 100644 index 02385f2..0000000 --- a/src/lasp/c/lasp_math_raw.c +++ /dev/null @@ -1,152 +0,0 @@ -// lasp_math_raw.c -// -// last-edit-by: J.A. de Jong -// -// Description: -// Operations working on raw arrays of floating point numbers -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-5) -#include "lasp_math_raw.h" -#if LASP_USE_BLAS == 1 -#include -#endif - -#ifdef __MKL_CBLAS_H__ -/* Symbol not present in MKL blas */ -#define blasint CBLAS_INDEX -#else -#endif - -void d_elem_prod_d(d res[], - const d arr1[], - const d arr2[], - const us size) { - - #if LASP_USE_BLAS == 1 - - #if LASP_DEBUG - - if(arr1 == arr2) { - DBGWARN("d_elem_prod_d: Array 1 and array 2 point to the same" - " memory. This results in pointer aliasing, for which" - " testing is still to be done. Results might be" - " unrealiable."); - } - - #endif - - - #if LASP_DOUBLE_PRECISION - #define elem_prod_fun cblas_dsbmv - #else - #define elem_prod_fun cblas_ssbmv - #endif - /* These parameters do not matter for this specific case */ - const CBLAS_ORDER mat_order= CblasColMajor; - const CBLAS_UPLO uplo = CblasLower; - - /* Extra multiplication factor */ - const d alpha = 1.0; - - /* void cblas_dsbmv(OPENBLAS_CONST enum CBLAS_ORDER order, */ - /* OPENBLAS_CONST enum CBLAS_UPLO Uplo, */ - /* OPENBLAS_CONST blasint N, */ - /* OPENBLAS_CONST blasint K, */ - /* OPENBLAS_CONST double alpha, */ - /* OPENBLAS_CONST double *A, */ - /* OPENBLAS_CONST blasint lda, */ - /* OPENBLAS_CONST double *X, */ - /* OPENBLAS_CONST blasint incX, */ - /* OPENBLAS_CONST double beta, */ - /* double *Y, */ - /* OPENBLAS_CONST blasint incY); */ - - elem_prod_fun(mat_order, - uplo, - (blasint) size, - 0, // Just the diagonal; 0 super-diagonal bands - alpha, /* Multiplication factor alpha */ - arr1, - 1, /* LDA */ - arr2, /* x */ - 1, /* incX = 1 */ - 0.0, /* Beta */ - res, /* The Y matrix to write to */ - 1); /* incY */ - #undef elem_prod_fun - - #else /* No blas routines, routine is very simple, but here we - * go! */ - DBGWARN("Performing slow non-blas vector-vector multiplication"); - for(us i=0;i -#include - -#if LASP_USE_BLAS == 1 -#include - -#elif LASP_USE_BLAS == 0 -#else -#error "LASP_USE_BLAS should be set to either 0 or 1" -#endif - - -/// Positive infinite -#define d_inf (INFINITY) - -#ifdef M_PI -static const d number_pi = M_PI; -#else -static const d number_pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679; -#endif - -/** - * Set all elements in an array equal to val - * - * @param to - * @param val - * @param size - */ -static inline void d_set(d to[],d val,us size) { - for(us i=0;ib?a:b; -} - - -/** - * Return the dot product of two arrays, one of them complex-valued, - * the other real-valued - * - * @param a the complex-valued array - * @param b the real-valued array - * @param size the size of the arrays. *Should be equal-sized!* - * - * @return the dot product - */ -static inline c cd_dot(const c a[],const d b[],us size){ - c result = 0; - us i; - for(i=0;i max) max=a[i]; - } - return max; -} -/** - * Compute the minimum of an array - * - * @param a array - * @param size size of the array - * @return minimum - */ -static inline d d_min(const d a[],us size){ - us i; - d min = a[0]; - for(i=1;i min) min=a[i]; - } - return min; -} - -/** - * Compute the \f$ L_2 \f$ norm of an array of doubles - * - * @param a Array - * @param size Size of array - */ -static inline d d_norm(const d a[],us size){ -#if LASP_USE_BLAS == 1 - return cblas_dnrm2(size,a,1); -#else - d norm = 0; - us i; - for(i=0;i - -/* #ifdef linux */ -#define _GNU_SOURCE -#include -#include -#ifndef MS_WIN64 -#include -#endif -/* #endif */ - -typedef struct { - void* job_ptr; - bool running; - bool ready; -} Job; - -typedef struct JobQueue_s { -#ifdef LASP_PARALLEL - pthread_mutex_t mutex; - pthread_cond_t cv_plus; /**< Condition variable for the - * "workers". */ - pthread_cond_t cv_minus; /**< Condition variable for the - * main thread. */ -#endif - Job* jobs; /**< Pointer to job vector */ - us max_jobs; /**< Stores the maximum number of - * items */ -} JobQueue; - -static us count_jobs(JobQueue* jq) { - fsTRACE(15); - us njobs = 0; - for(us i=0;imax_jobs;i++){ - if(jq->jobs[i].ready) - njobs++; - } - return njobs; -} -static Job* get_ready_job(JobQueue* jq) { - fsTRACE(15); - Job* j = jq->jobs; - for(us i=0;imax_jobs;i++){ - if(j->ready && !j->running) - return j; - j++; - } - return NULL; -} -void print_job_queue(JobQueue* jq) { - fsTRACE(15); - for(us i=0;imax_jobs;i++) { - printf("Job %zu", i); - if(jq->jobs[i].ready) - printf(" available"); - if(jq->jobs[i].running) - printf(" running"); - - printf(" - ptr %zu\n", (us) jq->jobs[i].job_ptr); - - } - feTRACE(15); -} - -#ifdef LASP_PARALLEL -#define LOCK_MUTEX \ - /* Lock the mutex to let the threads wait initially */ \ - int rv = pthread_mutex_lock(&jq->mutex); \ - if(rv !=0) { \ - WARN("Mutex lock failed"); \ - } -#define UNLOCK_MUTEX \ - rv = pthread_mutex_unlock(&jq->mutex); \ - if(rv !=0) { \ - WARN("Mutex unlock failed"); \ - } - -#else -#define LOCK_MUTEX -#define UNLOCK_MUTEX -#endif // LASP_PARALLEL - -JobQueue* JobQueue_alloc(const us max_jobs) { - TRACE(15,"JobQueue_alloc"); - /* if(max_jobs > LASP_MAX_NUM_CHANNELS) { */ - /* WARN("Max jobs restricted to LASP_MAX_NUM_CHANNELS"); */ - /* return NULL; */ - /* } */ - JobQueue* jq = a_malloc(sizeof(JobQueue)); - - - if(!jq) { - WARN("Allocation of JobQueue failed"); - return NULL; - } - jq->max_jobs = max_jobs; - - jq->jobs = a_malloc(max_jobs*sizeof(Job)); - if(!jq->jobs) { - WARN("Allocation of JobQueue jobs failed"); - return NULL; - } - - Job* j = jq->jobs; - for(us jindex=0;jindexjob_ptr = NULL; - j->ready = false; - j->running = false; - j++; - } - -#ifdef LASP_PARALLEL - /* Initialize thread mutex */ - int rv = pthread_mutex_init(&jq->mutex,NULL); - if(rv !=0) { - WARN("Mutex initialization failed"); - return NULL; - } - rv = pthread_cond_init(&jq->cv_plus,NULL); - if(rv !=0) { - WARN("Condition variable initialization failed"); - return NULL; - } - - rv = pthread_cond_init(&jq->cv_minus,NULL); - if(rv !=0) { - WARN("Condition variable initialization failed"); - return NULL; - } - -#endif // LASP_PARALLEL - /* print_job_queue(jq); */ - return jq; -} - -void JobQueue_free(JobQueue* jq) { - - TRACE(15,"JobQueue_free"); - dbgassert(jq,NULLPTRDEREF "jq in JobQueue_free"); - - int rv; - - if(count_jobs(jq) != 0) { - WARN("Job queue not empty!"); - } - - a_free(jq->jobs); - -#ifdef LASP_PARALLEL - /* Destroy the mutexes and condition variables */ - rv = pthread_mutex_destroy(&jq->mutex); - if(rv != 0){ - WARN("Mutex destroy failed. Do not know what to do."); - } - - rv = pthread_cond_destroy(&jq->cv_plus); - if(rv != 0){ - WARN("Condition variable destruction failed. " - "Do not know what to do."); - } - - rv = pthread_cond_destroy(&jq->cv_minus); - if(rv != 0){ - WARN("Condition variable destruction failed. " - "Do not know what to do."); - } -#endif // LASP_PARALLEL - -} - -int JobQueue_push(JobQueue* jq,void* job_ptr) { - - TRACE(15,"JobQueue_push"); - dbgassert(jq,NULLPTRDEREF "jq in JobQueue_push"); - - /* print_job_queue(jq); */ - /* uVARTRACE(15,(us) job_ptr); */ - - LOCK_MUTEX; - - us max_jobs = jq->max_jobs; - -#ifdef LASP_PARALLEL - /* Check if queue is full */ - while(count_jobs(jq) == max_jobs) { - - WARN("Queue full. Wait until some jobs are done."); - rv = pthread_cond_wait(&jq->cv_minus,&jq->mutex); - if(rv !=0) { - WARN("Condition variable wait failed"); - } - } -#else - /* If job queue is full, not in parallel, we just fail to add something - * without waiting*/ - if(count_jobs(jq) == max_jobs) { - return LASP_FAILURE; - } -#endif // LASP_PARALLEL - - dbgassert(count_jobs(jq) != max_jobs, - "Queue cannot be full!"); - - /* Queue is not full try to find a place, fill it */ - Job* j = jq->jobs; - us i; - for(i=0;iready == false ) { - dbgassert(j->job_ptr==NULL,"Job ptr should be 0"); - dbgassert(j->ready==false,"Job cannot be assigned"); - break; - } - j++; - } - dbgassert(i!=jq->max_jobs,"Should have found a job!"); - - j->job_ptr = job_ptr; - j->ready = true; - -#ifdef LASP_PARALLEL - /* Notify worker threads that a new job has arrived */ - if(count_jobs(jq) == max_jobs) { - /* Notify ALL threads. Action required! */ - rv = pthread_cond_broadcast(&jq->cv_plus); - if(rv !=0) { - WARN("Condition variable broadcast failed"); - } - - } else { - /* Notify some thread that there has been some change to - * the Queue */ - rv = pthread_cond_signal(&jq->cv_plus); - if(rv !=0) { - WARN("Condition variable signal failed"); - } - } -#endif // LASP_PARALLEL - - /* print_job_queue(jq); */ - - UNLOCK_MUTEX; - - return LASP_SUCCESS; -} -void* JobQueue_assign(JobQueue* jq) { - - TRACE(15,"JobQueue_assign"); - - LOCK_MUTEX; - - Job* j; -#ifdef LASP_PARALLEL - /* Wait until a job is available */ - while ((j=get_ready_job(jq))==NULL) { - - TRACE(15,"JobQueue_assign: no ready job"); - pthread_cond_wait(&jq->cv_plus,&jq->mutex); - - } -#else - if(count_jobs(jq) == 0) { return NULL; } - else { j = get_ready_job(jq); } -#endif // LASP_PARALLEL - - TRACE(16,"JobQueue_assign: found ready job. Assigned to:"); -#ifdef LASP_DEBUG -#ifdef LASP_PARALLEL - pthread_t thisthread = pthread_self(); - iVARTRACE(16,thisthread); -#endif -#endif - - /* print_job_queue(jq); */ - /* Find a job from the queue, assign it and return it */ - j->running = true; - -#ifdef LASP_PARALLEL - if(count_jobs(jq) > 1) { - /* Signal different thread that there is more work to do */ - rv = pthread_cond_signal(&jq->cv_plus); - if(rv !=0) { - WARN("Condition variable broadcast failed"); - } - } -#endif - - UNLOCK_MUTEX; - - TRACE(15,"End JobQueue_assign"); - - return j->job_ptr; -} -void JobQueue_done(JobQueue* jq,void* job_ptr) { - - TRACE(15,"JobQueue_done"); - dbgassert(jq,NULLPTRDEREF "jq in JobQueue_done"); - - LOCK_MUTEX; - - /* print_job_queue(jq); */ - - /* Find the job from the queue, belonging to the job_ptr */ - Job* j=jq->jobs; - us i; - for(i=0;imax_jobs;i++) { - iVARTRACE(10,i); - if(j->ready && j->running && j->job_ptr == job_ptr) { - TRACE(15,"Found the job that has been done:"); - j->ready = false; - j->job_ptr = NULL; - j->running = false; - break; - } - j++; - } - - /* print_job_queue(jq); */ - -#ifdef LASP_PARALLEL - /* Job done, broadcast this */ - rv = pthread_cond_signal(&jq->cv_minus); - if(rv !=0) { - WARN("Condition variable broadcast failed"); - } -#endif - - UNLOCK_MUTEX; -} - -#ifdef LASP_PARALLEL -void JobQueue_wait_alldone(JobQueue* jq) { - TRACE(15,"JobQueue_wait_alldone"); - dbgassert(jq,NULLPTRDEREF "jq in JobQueue_wait_alldone"); - - LOCK_MUTEX; - - /* Wait until number of jobs is 0 */ - while (count_jobs(jq)!=0) { - - if(rv !=0) { - WARN("Condition variable broadcast failed"); - } - - pthread_cond_wait(&jq->cv_minus,&jq->mutex); - } - - UNLOCK_MUTEX; - -} -#endif - - -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_mq.h b/src/lasp/c/lasp_mq.h deleted file mode 100644 index e183229..0000000 --- a/src/lasp/c/lasp_mq.h +++ /dev/null @@ -1,78 +0,0 @@ -// mq.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// Multithreaded job queue implementation -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef MQ_H -#define MQ_H -#include "lasp_types.h" - -typedef struct JobQueue_s JobQueue; - -/** - * Allocate a new job queue. - * - * @param max_msg Maximum number of jobs that can be put in the - * queue. - * - * @return Pointer to new JobQueue instance. NULL on error. - */ -JobQueue* JobQueue_alloc(const us max_msg); - -/** - * Free an existing job queue. If it is not empty and threads are - * still waiting for jobs, the behaviour is undefined. So please - * make sure all threads are done before free'ing the queue. - * - * @param jq: JobQueue to free - */ -void JobQueue_free(JobQueue* jq); - -/** - * Pops a job from the queue. Waits indefinitely until some job is - * available. If in parallel mode. In serial mode, it returns NULL if no job - * is available (LASP_PARALLEL compilation flag not set). - * - * @param jq: JobQueue handle - * @return Pointer to the job, NULL on error. - */ -void* JobQueue_assign(JobQueue* jq); - -/** - * Tell the queue the job that has been popped is done. Only after - * this function call, the job is really removed from the queue. - * - * @param jq: JobQueue handle - * @param job - */ -void JobQueue_done(JobQueue* jq,void* job); - -/** - * A push on the job queue will notify one a single thread that is - * blocked waiting in the JobQueue_assign() function. If the job - * queue is full, however all waiters will be signaled and the - * function will block until there is some space in the job queue. - * - * @param jp JobQueue - * @param job_ptr Pointer to job to be done - * @return 0 on success. - */ -int JobQueue_push(JobQueue* jq,void* job_ptr); - -/** - * Wait until the job queue is empty. Please use this function with - * caution, as it will block indefinitely in case the queue never gets - * empty. The purpose of this function is to let the main thread wait - * until all task workers are finished. - * - */ - -#ifdef LASP_PARALLEL -void JobQueue_wait_alldone(JobQueue*); -#endif // LASP_PARALLEL - -#endif // MQ_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_nprocs.c b/src/lasp/c/lasp_nprocs.c deleted file mode 100644 index 4467499..0000000 --- a/src/lasp/c/lasp_nprocs.c +++ /dev/null @@ -1,18 +0,0 @@ -#include "lasp_nprocs.h" -#ifdef MS_WIN64 -#include -#else -// Used for obtaining the number of processors -#include -#endif - -us getNumberOfProcs() { -#if MS_WIN64 -// https://stackoverflow.com/questions/150355/programmatically-find-the-number-of-cores-on-a-machine - SYSTEM_INFO sysinfo; - GetSystemInfo(&sysinfo); - return sysinfo.dwNumberOfProcessors; -#else // Linux, easy - return get_nprocs(); -#endif -} diff --git a/src/lasp/c/lasp_nprocs.h b/src/lasp/c/lasp_nprocs.h deleted file mode 100644 index ea2f34c..0000000 --- a/src/lasp/c/lasp_nprocs.h +++ /dev/null @@ -1,25 +0,0 @@ -// lasp_nprocs.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: Implemententation of a function to determine the number -// of processors. Useful for spreading out computational load over multiple -// cores. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_NPROCS_H -#define LASP_NPROCS_H -#include "lasp_types.h" -#ifdef __cplusplus -extern C { -#endif - - /** - * @return The number of SMP processors - */ - us getNumberOfProcs(); - -#ifdef __cplusplus -} -#endif -#endif // LASP_NPROCS_H diff --git a/src/lasp/c/lasp_ps.c b/src/lasp/c/lasp_ps.c deleted file mode 100644 index c8e664d..0000000 --- a/src/lasp/c/lasp_ps.c +++ /dev/null @@ -1,181 +0,0 @@ -// lasp_ps.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-5) -#include "lasp_ps.h" -#include "lasp_fft.h" -#include "lasp_alloc.h" -#include "lasp_alg.h" -#include "lasp_assert.h" - -typedef struct PowerSpectra_s { - - vd window; - d win_pow; /**< The power of the window */ - Fft* fft; /**< Fft routines storage */ -} PowerSpectra; - -PowerSpectra* PowerSpectra_alloc(const us nfft, - const WindowType wt) { - - fsTRACE(15); - int rv; - - /* Check nfft */ - if(nfft % 2 != 0) { - WARN("nfft should be even"); - return NULL; - } - - /* ALlocate space */ - Fft* fft = Fft_create(nfft); - if(fft == NULL) { - WARN("Fft allocation failed"); - return NULL; - } - - PowerSpectra* ps = a_malloc(sizeof(PowerSpectra)); - if(!ps) { - WARN("Allocation of PowerSpectra memory failed"); - Fft_free(fft); - return NULL; - } - ps->fft = fft; - - /* Allocate vectors and matrices */ - ps->window = vd_alloc(nfft); - - rv = window_create(wt,&ps->window,&ps->win_pow); - check_overflow_vx(ps->window); - if(rv!=0) { - WARN("Error creating window function, continuing anyway"); - } - feTRACE(15); - return ps; -} - -void PowerSpectra_free(PowerSpectra* ps) { - fsTRACE(15); - Fft_free(ps->fft); - vd_free(&ps->window); - a_free(ps); - feTRACE(15); -} - - -void PowerSpectra_compute(const PowerSpectra* ps, - const dmat * timedata, - cmat * result) { - - fsTRACE(15); - - dbgassert(ps && timedata && result,NULLPTRDEREF); - - const us nchannels = timedata->n_cols; - const us nfft = Fft_nfft(ps->fft); - uVARTRACE(15,nchannels); - const d win_pow = ps->win_pow; - dVARTRACE(15,win_pow); - - /* Sanity checks for the matrices */ - dbgassert(timedata->n_rows == nfft,"timedata n_rows " - "should be equal to nfft"); - - dbgassert(result->n_rows == nfft/2+1,"result n_rows " - "should be equal to nfft/2+1"); - - dbgassert(result->n_cols == nchannels*nchannels,"result n_cols " - "should be equal to nchannels*nchannels"); - - - /* Multiply time data with the window and store result in - * timedata_work. */ - dmat timedata_work = dmat_alloc(nfft,nchannels); - for(us i=0;iwindow); - - vd_free(&column); - vd_free(&column_work); - } - check_overflow_xmat(timedata_work); - /* print_dmat(&timedata_work); */ - - /* Compute fft of the time data */ - cmat fft_work = cmat_alloc(nfft/2+1,nchannels); - Fft_fft(ps->fft, - &timedata_work, - &fft_work); - - dmat_free(&timedata_work); - - TRACE(15,"fft done"); - - /* Scale fft such that power is easily comxputed */ - const c scale_fac = d_sqrt(2/win_pow)/nfft; - - /* POWER SPECTRAL DENSITY? Scale fft such that power is easily computed - * - Multiply power spectral density by 2 except at f=0 and f=fNq - * - Divide by energy of window function = nfft * window_power - * - .. sqrt(factors) because it is applied to output fft instead of psd */ - // const c scale_fac = d_sqrt(2/(nfft*win_pow)); - - cmat_scale(&fft_work,scale_fac); - TRACE(15,"scale done"); - - for(us i=0;i< nchannels;i++) { - /* Multiply DC term by 1/sqrt(2) */ - *getcmatval(&fft_work,0,i) *= 1/d_sqrt(2.)+0*I; - - /* Multiply Nyquist term by 1/sqrt(2) */ - *getcmatval(&fft_work,nfft/2,i) *= 1/d_sqrt(2.)+0*I; - } - check_overflow_xmat(fft_work); - - /* print_cmat(&fft_work); */ - TRACE(15,"Nyquist and DC correction done"); - - vc j_vec_conj = vc_alloc(nfft/2+1); - - /* Compute Cross-power spectra and store result */ - for(us i =0; i -#include - -#if LASP_DOUBLE_PRECISION == 1 -#define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT64 -#define LASP_NUMPY_COMPLEX_TYPE NPY_COMPLEX128 -#else -#define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT32 -#endif -#include - -/** - * Function passed to Python to use for cleanup of - * foreignly obtained data. - **/ -#define LASP_CAPSULE_NAME "pyarray_data_destructor" -static inline void capsule_cleanup(PyObject *capsule) { - void *memory = PyCapsule_GetPointer(capsule, LASP_CAPSULE_NAME); - free(memory); -} - -/** - * Create a numpy array from a raw data pointer. - * - * @param data pointer to data, assumes d* for data - * @param transfer_ownership If set to true, the created Numpy array will be - * responsible for freeing the data. - * @param F_contiguous It set to true, the data is assumed to be column-major - * ordered in memory (which means we can find element d[r,c] by d[n_cols*r+c], - * if set to false. Data is assumed to be row-major ordered (which means we - * find element d[r,c] by d[n_rows*c+r] - * - * @return Numpy array - */ -static inline PyObject *data_to_ndarray(void *data, int n_rows, int n_cols, - int typenum, bool transfer_ownership, - bool F_contiguous) { - - /* fprintf(stderr, "Enter data_to_ndarray\n"); */ - assert(data); - - PyArray_Descr *descr = PyArray_DescrFromType(typenum); - if (!descr) - return NULL; - - npy_intp dims[2] = {n_rows, n_cols}; - npy_intp strides[2]; - - int flags = 0; - if (F_contiguous) { - flags |= NPY_ARRAY_FARRAY; - strides[0] = descr->elsize; - strides[1] = descr->elsize * n_rows; - } else { - strides[0] = descr->elsize * n_rows; - strides[1] = descr->elsize; - } - - PyArrayObject *arr = - (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, - descr, // Description - 2, // nd - dims, // dimensions - strides, // strides - data, // Data pointer - flags, // Flags - NULL // obj - ); - - if (!arr) { - fprintf(stderr, "arr = 0!"); - return NULL; - } - - if (transfer_ownership) { - // The default destructor of Python cannot free the data, as it is allocated - // with malloc. Therefore, with this code, we tell Numpy/Python to use - // the capsule_cleanup constructor. See: - // https://stackoverflow.com/questions/54269956/crash-of-jupyter-due-to-the-use-of-pyarray-enableflags/54278170#54278170 - // Note that in general it was disadvised to build all C code with MinGW on - // Windows. We do it anyway, see if we find any problems on the way. - /* PyObject *capsule = PyCapsule_New(data, LASP_CAPSULE_NAME, capsule_cleanup); */ - /* int res = PyArray_SetBaseObject(arr, capsule); */ - /* if (res != 0) { */ - /* fprintf(stderr, "Failed to set base object of array!"); */ - /* return NULL; */ - /* } */ - - PyArray_ENABLEFLAGS(arr, NPY_OWNDATA); - } - - return (PyObject *)arr; -} - -#undef LASP_CAPSULE_NAME -#endif // LASP_PYARRAY_H diff --git a/src/lasp/c/lasp_python.h b/src/lasp/c/lasp_python.h deleted file mode 100644 index 27dc8d9..0000000 --- a/src/lasp/c/lasp_python.h +++ /dev/null @@ -1,59 +0,0 @@ -// ascee_python.h -// -// Author: J.A. de Jong - ASCEE - Redu-Sone -// -// Description: -// Some routines to generate numpy arrays from matrices and vectors. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_PYTHON_H -#define LASP_PYTHON_H - -#ifdef __cplusplus -#error "Cannot compile this file with C++" -#endif -#include "lasp_types.h" -#include "lasp_mat.h" -#include "lasp_pyarray.h" - -/** - * @brief Create a Numpy array from a dmat structure - * - * @param mat Pointer to the dmat - * @param transfer_ownership If set to true, the created Numpy array will - * obtain ownership of the data and calls free() on destruction. The dmat will - * have the foreign_data flag set, such that it won't free the data. - * - * @return - */ -static inline PyObject* dmat_to_ndarray(dmat* mat,bool transfer_ownership) { - dbgassert(mat,NULLPTRDEREF); - dbgassert(mat->_data,NULLPTRDEREF); - /* fprintf(stderr, "Enter dmat_to_ndarray\n"); */ - - // Dimensions given in wrong order, as mat is - // Fortran-contiguous. Later on we transpose the result. This is - // more easy than using the PyArray_New syntax. - PyObject* arr = data_to_ndarray(mat->_data, - mat->n_rows, - mat->n_cols, - LASP_NUMPY_FLOAT_TYPE, - transfer_ownership, - true); // Fortran-contiguous - - if(transfer_ownership) { - mat->_foreign_data = true; - } - - if(!arr) { - WARN("Array creation failure"); - feTRACE(15); - return NULL; - } - /* fprintf(stderr, "Exit dmat_to_ndarray\n"); */ - - return arr; -} - -#endif // LASP_PYTHON_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_signals.h b/src/lasp/c/lasp_signals.h deleted file mode 100644 index 71205f0..0000000 --- a/src/lasp/c/lasp_signals.h +++ /dev/null @@ -1,32 +0,0 @@ -// lasp_signals.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// Several signal functions -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_SIGNALS_H -#define LASP_SIGNALS_H -#include "lasp_mat.h" - -/** - * Compute the signal power, that is \f$ \frac{1}{N} \sum_{i=0}^{N-1} - * v_i^2 \f$ - * - * @param[in] signal Signal to compute the power of. - * @return the signal power - */ -static inline d signal_power(vd* signal) { - d res = 0; - for(us i=0;in_rows;i++) { - res+= d_pow(*getvdval(signal,i),2); - } - res /= signal->n_rows; - return res; -} - - - -#endif // LASP_SIGNALS_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_slm.c b/src/lasp/c/lasp_slm.c deleted file mode 100644 index ccbdb82..0000000 --- a/src/lasp/c/lasp_slm.c +++ /dev/null @@ -1,293 +0,0 @@ -#define TRACERPLUS (-5) -#include "lasp_slm.h" -#include "lasp_assert.h" -#include "lasp_tracer.h" - -typedef struct Slm { - Sosfilterbank *prefilter; /// Pre-filter, A, or C. If NULL, not used. - Sosfilterbank *bandpass; /// Filterbank. If NULL, not used - Sosfilterbank **splowpass; /// Used for time-weighting of the squared signal - d ref_level; /// Reference value for computing decibels - us downsampling_fac; /// Every x'th sample is returned. - us cur_offset; /// Storage for offset point in input arrays - vd Pm; /// Storage for the computing the mean of the square of the signal. - vd Pmax; /// Storage for maximum computed signal power so far. - vd Ppeak; /// Storage for computing peak powers so far. - us N; /// Counter for the number of time samples counted that came in - -} Slm; - -Slm *Slm_create(Sosfilterbank *prefilter, Sosfilterbank *bandpass, const d fs, - const d tau, const d ref_level, us *downsampling_fac) { - fsTRACE(15); - assertvalidptr(downsampling_fac); - - Slm *slm = NULL; - if (ref_level <= 0) { - WARN("Invalid reference level"); - return NULL; - } else if (fs <= 0) { - WARN("Invalid sampling frequency"); - return NULL; - } - - slm = (Slm *)a_malloc(sizeof(Slm)); - slm->ref_level = ref_level; - slm->prefilter = prefilter; - slm->bandpass = bandpass; - - /// Compute the downsampling factor. This one is chosen based on the - /// lowpass filter. Which has a -3 dB point of f = 1/(tau*2*pi). See LASP - /// documentation for the computation of its minus 20 dB point. We set the - /// reduction in its 'sampling frequency' such that its noise is at a level - /// of 20 dB less than its 'signal'. - us ds_fac; - if (tau > 0) { - // A reasonable 'framerate' for the sound level meter, based on the - // filtering time constant. - d fs_slm = 10 / tau; - dVARTRACE(15, fs_slm); - if(fs_slm < 30) { - fs_slm = 30; - } - ds_fac = (us)(fs / fs_slm); - if (ds_fac == 0) { - // If we get 0, it should be 1 - ds_fac++; - } - } else { - ds_fac = 1; - } - slm->downsampling_fac = ds_fac; - *downsampling_fac = ds_fac; - slm->cur_offset = 0; - - /// Create the single pole lowpass - us filterbank_size; - if (bandpass) { - filterbank_size = Sosfilterbank_getFilterbankSize(bandpass); - } else { - filterbank_size = 1; - } - - if (tau > 0) { - - vd lowpass_sos = vd_alloc(6); - d b0 = 1.0 / (1 + 2 * tau * fs); - *getvdval(&lowpass_sos, 0) = b0; - *getvdval(&lowpass_sos, 1) = b0; - *getvdval(&lowpass_sos, 2) = 0; - *getvdval(&lowpass_sos, 3) = 1; - *getvdval(&lowpass_sos, 4) = (1 - 2 * tau * fs) * b0; - *getvdval(&lowpass_sos, 5) = 0; - - slm->splowpass = a_malloc(filterbank_size * sizeof(Sosfilterbank *)); - for (us ch = 0; ch < filterbank_size; ch++) { - /// Allocate a filterbank with one channel and one section. - slm->splowpass[ch] = Sosfilterbank_create(0, 1, 1); - Sosfilterbank_setFilter(slm->splowpass[ch], 0, lowpass_sos); - } - vd_free(&lowpass_sos); - } else { - /// No low-pass filtering. Tau set to zero - slm->splowpass = NULL; - } - - /// Initialize statistics gatherers - slm->Ppeak = vd_alloc(filterbank_size); - slm->Pmax = vd_alloc(filterbank_size); - slm->Pm = vd_alloc(filterbank_size); - slm->N = 0; - vd_set(&(slm->Ppeak), 0); - vd_set(&(slm->Pmax), 0); - vd_set(&(slm->Pm), 0); - feTRACE(15); - return slm; -} - -dmat Slm_run(Slm *slm, vd *input_data) { - fsTRACE(15); - assertvalidptr(slm); - assert_vx(input_data); - - /// First step: run the input data through the pre-filter - vd prefiltered; - if (slm->prefilter) - prefiltered = Sosfilterbank_filter(slm->prefilter, input_data); - else { - prefiltered = dmat_foreign(input_data); - } - dmat bandpassed; - if (slm->bandpass) { - bandpassed = Sosfilterbank_filter(slm->bandpass, &prefiltered); - } else { - bandpassed = dmat_foreign(&prefiltered); - } - us filterbank_size = bandpassed.n_cols; - - /// Next step: square all values. We do this in-place. Then we filter for - /// each channel. - d ref_level = slm->ref_level; - d *tmp; - - /// Pre-calculate the size of the output data - us downsampling_fac = slm->downsampling_fac; - us samples_bandpassed = bandpassed.n_rows; - iVARTRACE(15, samples_bandpassed); - iVARTRACE(15, downsampling_fac); - us cur_offset = slm->cur_offset; - - /// Compute the number of samples output - us nsamples_output = samples_bandpassed; - if (downsampling_fac > 1) { - nsamples_output = (samples_bandpassed - cur_offset) / downsampling_fac; - if(nsamples_output > samples_bandpassed) { - // This means overflow of unsigned number calculations - nsamples_output = 0; - } - while(nsamples_output * downsampling_fac + cur_offset < samples_bandpassed) { - nsamples_output++; - } - } - - iVARTRACE(15, nsamples_output); - iVARTRACE(15, cur_offset); - dmat levels = dmat_alloc(nsamples_output, filterbank_size); - us N, ch; - - for (ch = 0; ch < bandpassed.n_cols; ch++) { - iVARTRACE(15, ch); - vd chan = dmat_column(&bandpassed, ch); - /// Inplace squaring of the signal - for (us sample = 0; sample < bandpassed.n_rows; sample++) { - tmp = getdmatval(&bandpassed, sample, ch); - *tmp = *tmp * *tmp; - *getvdval(&(slm->Ppeak), ch) = d_max(*getvdval(&(slm->Ppeak), ch), *tmp); - } - - // Now that all data for the channel is squared, we can run it through - // the low-pass filter - cur_offset = slm->cur_offset; - - /// Apply single-pole lowpass filter for current filterbank channel - TRACE(15, "Start filtering"); - vd power_filtered; - if (slm->splowpass) { - power_filtered = Sosfilterbank_filter(slm->splowpass[ch], &chan); - } else { - power_filtered = dmat_foreign(&chan); - } - TRACE(15, "Filtering done"); - dbgassert(chan.n_rows == power_filtered.n_rows, "BUG"); - - /// Output resulting levels at a lower interval - us i = 0; - N = slm->N; - d *Pm = getvdval(&(slm->Pm), ch); - while (cur_offset < samples_bandpassed) { - iVARTRACE(10, i); - iVARTRACE(10, cur_offset); - - /// Filtered power. - const d P = *getvdval(&power_filtered, cur_offset); - dVARTRACE(15, P); - - /// Compute maximum, compare to current maximum - *getvdval(&(slm->Pmax), ch) = d_max(*getvdval(&(slm->Pmax), ch), P); - - /// Update mean power - d Nd = (d) N; - *Pm = (*Pm*Nd + P ) / (Nd+1); - N++; - dVARTRACE(15, *Pm); - - /// Compute level - d level = 10 * d_log10((P + d_epsilon ) / ref_level / ref_level); - - *getdmatval(&levels, i++, ch) = level; - cur_offset = cur_offset + downsampling_fac; - } - iVARTRACE(15, cur_offset); - iVARTRACE(15, i); - dbgassert(i == nsamples_output, "BUG"); - - vd_free(&power_filtered); - vd_free(&chan); - } - /// Update sample counter - dbgassert(ch >0, "BUG"); - slm->N = N; - slm->cur_offset = cur_offset - samples_bandpassed; - - vd_free(&prefiltered); - dmat_free(&bandpassed); - feTRACE(15); - return levels; -} - - -static inline vd levels_from_power(const vd* power,const d ref_level){ - fsTRACE(15); - - vd levels = dmat_alloc_from_dmat(power); - for(us i=0; i< levels.n_rows; i++) { - *getvdval(&levels, i) = 10 * d_log10( - (*getvdval(power, i) + d_epsilon) / ref_level / ref_level); - - } - feTRACE(15); - return levels; -} - -vd Slm_Lpeak(Slm* slm) { - fsTRACE(15); - assertvalidptr(slm); - vd Lpeak = levels_from_power(&(slm->Ppeak), slm->ref_level); - feTRACE(15); - return Lpeak; -} - -vd Slm_Lmax(Slm* slm) { - fsTRACE(15); - assertvalidptr(slm); - vd Lmax = levels_from_power(&(slm->Pmax), slm->ref_level); - feTRACE(15); - return Lmax; -} - -vd Slm_Leq(Slm* slm) { - fsTRACE(15); - assertvalidptr(slm); - print_vd(&(slm->Pm)); - vd Leq = levels_from_power(&(slm->Pm), slm->ref_level); - feTRACE(15); - return Leq; -} - -void Slm_free(Slm *slm) { - fsTRACE(15); - assertvalidptr(slm); - if (slm->prefilter) { - Sosfilterbank_free(slm->prefilter); - } - - us filterbank_size; - if (slm->bandpass) { - filterbank_size = Sosfilterbank_getFilterbankSize(slm->bandpass); - Sosfilterbank_free(slm->bandpass); - } else { - filterbank_size = 1; - } - if (slm->splowpass) { - for (us ch = 0; ch < filterbank_size; ch++) { - Sosfilterbank_free(slm->splowpass[ch]); - } - a_free(slm->splowpass); - } - vd_free(&(slm->Ppeak)); - vd_free(&(slm->Pmax)); - vd_free(&(slm->Pm)); - a_free(slm); - - feTRACE(15); -} diff --git a/src/lasp/c/lasp_slm.h b/src/lasp/c/lasp_slm.h deleted file mode 100644 index 479b383..0000000 --- a/src/lasp/c/lasp_slm.h +++ /dev/null @@ -1,96 +0,0 @@ -// lasp_slm.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: Multi-purpose implementation of a real time Sound Level Meter, -// can be used for full-signal filtering, (fractional) octave band filtering, -// etc. -// ////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_SLM_H -#define LASP_SLM_H -#include "lasp_types.h" -#include "lasp_mat.h" -#include "lasp_sosfilterbank.h" - -typedef struct Slm Slm; - -#define TAU_FAST (1.0/8.0) -#define TAU_SLOW (1.0) -#define TAU_IMPULSE (35e-3) - -/** - * Initializes a Sound level meter. NOTE: Sound level meter takes over - * ownership of pointers to the filterbanks for prefiltering and band-pass - * filtering! After passing them to the constructor of the Slm, they should not - * be touched anymore. - * - * @param[in] weighting: Sosfiterbank handle, used to pre-filter data. This is - * in most cases an A-weighting filter, or C-weighting. If NULL, no - * pre-filtering is done. That can be the case in situations of Z-weighting. - * @param[in] fb: Sosfiterbank handle, bandpass filters. - * @param[in] fs: Sampling frequency in [Hz], used for computing the - * downsampling factor and size of output arrays. - * @param[in] tau: Time constant of single pole low pass filter for squared data. - * Three pre-defined values can be used: TAU_FAST, for fast filtering, - * TAU_SLOW, for slow filtering, TAU_IMPULSE for impulse filtering - * @param[in] ref_level: Reference level when computing dB's. I.e. P_REF_AIR for - * sound pressure levels in air - * @param[out] downsampling_fac: Here, the used downsampling factor is stored - * which is used on returning data. If the value is for example 10, the - * 'sampling' frequency of output data from `Slm_run` is 4800 is fs is set to - * 48000. This downsampling factor is a function of the used time weighting. - * - * @return Slm: Handle of the sound level meter, NULL on error. - */ -Slm* Slm_create(Sosfilterbank* weighting,Sosfilterbank* bandpass, - const d fs, - const d tau, - const d ref_level, - us* downsampling_fac); - -/** - * Run the sound level meter on a piece of time data. - * - * @param[in] slm: Slm handle - * @param[in] input_data: Vector of input data samples. - * - * @return Output result of Sound Level values in [dB], for each bank in the filter - * bank. - */ -dmat Slm_run(Slm* slm, - vd* input_data); -/** - * Cleans up an existing Slm - * - * @param f slm handle - */ -void Slm_free(Slm* f); - -/** - * Returns the (raw) peak level computed so far. - * - * @param[in] slm: Slm handle - * @return Vector of peak level for each channel in the filterbank. - */ -vd Slm_Lpeak(Slm* slm); - -/** - * Returns the equivalent level computed so far. - * - * @param[in] slm: Slm handle - * @return Vector of equivalent levels for each channel in the filterbank. - */ -vd Slm_Leq(Slm* slm); - -/** - * Returns the maximum level computed so far. - * - * @param[in] slm: Slm handle - * @return Vector of maximum levels for each channel in the filterbank. - */ -vd Slm_Lmax(Slm* slm); - - -#endif // LASP_SLM_H -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_sosfilterbank.c b/src/lasp/c/lasp_sosfilterbank.c deleted file mode 100644 index 7b42147..0000000 --- a/src/lasp/c/lasp_sosfilterbank.c +++ /dev/null @@ -1,264 +0,0 @@ -#define TRACERPLUS (-5) -#include "lasp_sosfilterbank.h" -#include "lasp_mq.h" -#include "lasp_worker.h" -#include "lasp_nprocs.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; -#ifdef LASP_PARALLEL - JobQueue* jq; - Workers* workers; -#endif // LASP_PARALLEL - - -} Sosfilterbank; - -us Sosfilterbank_getFilterbankSize(const Sosfilterbank* fb) { - fsTRACE(15); - assertvalidptr(fb); - return fb->filterbank_size; - feTRACE(15); - -} - -int filter_single(void* worker_data,void* job); - -static inline us min(us a, us b){ - return ab?a:b; -} - -Sosfilterbank* Sosfilterbank_create( - const us nthreads_, - 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); - -#ifdef LASP_PARALLEL - fb->jq = NULL; - fb->workers = NULL; - dbgassert(nthreads_ <= LASP_MAX_NUM_THREADS, "Illegal number of threads"); - us nthreads; - us nprocs = getNumberOfProcs(); - - if(nthreads_ == 0) { - nthreads = min(max(nprocs/2,1), filterbank_size); - } else { - nthreads = nthreads_; - } - iVARTRACE(15, nthreads); - - if(nthreads > 1) { - if(!(fb->jq = JobQueue_alloc(filterbank_size))) { - Sosfilterbank_free(fb); - feTRACE(15); - return NULL; - } - - if(!(fb->workers = Workers_create(nthreads, - fb->jq, - NULL, - &filter_single, - NULL, - NULL))) { - Sosfilterbank_free(fb); - feTRACE(15); - return NULL; - } - } - -#endif // LASP_PARALLEL - 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); - iVARTRACE(15, filter_coefs.n_rows); - iVARTRACE(15, filter_no); - 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;indexsos)); - dmat_free(&(fb->state)); - -#ifdef LASP_PARALLEL - if(fb->workers) Workers_free(fb->workers); - if(fb->jq) JobQueue_free(fb->jq); -#endif // LASP_PARALLEL - - a_free(fb); - feTRACE(15); - -} - -typedef struct { - us filter_no; - us nsections; - us nsamples; - dmat sos; - dmat state; - dmat ys; -} Job; - -int filter_single(void* worker_data, void* job_) { - fsTRACE(15); - Job* job = (Job*) job_; - - us nsections = job->nsections; - us nsamples = job->nsamples; - dmat sos = job->sos; - - - for(us section=0;sectionstate),job->filter_no,section*2); - d w2 = *getdmatval(&(job->state),job->filter_no,section*2+1); - - d b0 = *getdmatval(&sos,job->filter_no,section*6+0); - d b1 = *getdmatval(&sos,job->filter_no,section*6+1); - d b2 = *getdmatval(&sos,job->filter_no,section*6+2); - /* d a0 = *getdmatval(&sos,job->filter_no,section*6+3); */ - d a1 = *getdmatval(&sos,job->filter_no,section*6+4); - d a2 = *getdmatval(&sos,job->filter_no,section*6+5); - - d* y = getdmatval(&(job->ys), 0, job->filter_no); - - for(us sample=0;samplestate),job->filter_no,section*2) = w1; - *getdmatval(&(job->state),job->filter_no,section*2+1) = w2; - - } - - feTRACE(15); - return 0; -} - - -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;filterworkers) { - assertvalidptr(fb->jq); - JobQueue_push(fb->jq, &(jobs[filter])); - } else { - -#endif // LASP_PARALLEL - /* No workers, we have to do it ourselves */ - filter_single(NULL,(void*) &(jobs[filter])); -#ifdef LASP_PARALLEL - } -#endif // LASP_PARALLEL - } -#ifdef LASP_PARALLEL - if(fb->workers) { - JobQueue_wait_alldone(fb->jq); - } -#endif // LASP_PARALLEL - feTRACE(15); - return ys; -} diff --git a/src/lasp/c/lasp_sosfilterbank.h b/src/lasp/c/lasp_sosfilterbank.h deleted file mode 100644 index 67c1815..0000000 --- a/src/lasp/c/lasp_sosfilterbank.h +++ /dev/null @@ -1,81 +0,0 @@ -// lasp_sosfilterbank.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: Implemententation of a discrete parallel filterbank using -// cascaded second order sections (sos) for each filter in the filterbank. In -// parallel mode, the filters are allocated over a set of Worker threads that -// actually perform the computation. In serial mode, all filters in the bank -// are computed in series. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_FILTERBANK_H -#define LASP_FILTERBANK_H -#include "lasp_config.h" -#include "lasp_types.h" -#include "lasp_mat.h" - -#define MAX_SOS_FILTER_BANK_SIZE 40 -#define MAX_SOS_FILTER_BANK_NSECTIONS 10 - -typedef struct Sosfilterbank Sosfilterbank; - -/** - * Initializes a Sosfilterbank. Sets all coefficients in such a way that the - * filter effectively does nothing (unit impulse response). - * @param[in] nthreads: The number of threads used in computation should be 1 - * <= nthreads <= LASP_MAX_NUM_THREADS. If set to 0, a sensible value is - * computed based on the number of filterbanks. - * @param[in] filterbank_size: The number of parallel filters in the bank - * - * @return Sosfilterbank handle - */ -Sosfilterbank* Sosfilterbank_create(const us nthreads, - const us filterbank_size, - const us nsections); - -/** - * Returns the number of channels in the filterbank (the filberbank size). - * - * @param[in] fb: Filterbank handle - * @return The number of filters in the bank - * */ -us Sosfilterbank_getFilterbankSize(const Sosfilterbank* fb); - -/** - * 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 b are the numerator coefficients and a 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 (and is equal to the number of rows in the output). - */ -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 -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_worker.c b/src/lasp/c/lasp_worker.c deleted file mode 100644 index 7a26a77..0000000 --- a/src/lasp/c/lasp_worker.c +++ /dev/null @@ -1,206 +0,0 @@ -// lasp_worker.c -// -// Author: J.A. de Jong -ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-5) - -/* The code in this file is not of any use for systems that are not able to do - * simultaneous multithreading, or should be adjusted for it. It is therefore - * only compiled in case LASP_PARALLEL flag is set */ -#include "lasp_worker.h" -#ifdef LASP_PARALLEL -#include "lasp_mq.h" -#include "lasp_alloc.h" -#include -#include "lasp_assert.h" -#include "lasp_tracer.h" - -typedef struct Workers_s { - JobQueue* jq; - worker_alloc_function w_alloc_fcn; - worker_function fn; - worker_free_function w_free_fcn; - - pthread_mutex_t global_data_mutex; - void* global_data; - - pthread_t worker_threads[LASP_MAX_NUM_THREADS]; - us num_workers; -} Workers; - -static void* threadfcn(void* data); - -Workers* Workers_create(const us num_workers, - JobQueue* jq, - worker_alloc_function init_fn, - worker_function fn, - worker_free_function free_fn, - void* thread_global_data) { - - TRACE(15,"Workers_create"); - - if(num_workers > LASP_MAX_NUM_THREADS) { - WARN("Number of workers too high in Workers_create"); - return NULL; - } - - /* dbgassert(init_fn,NULLPTRDEREF "init_fn"); */ - dbgassert(fn,NULLPTRDEREF "fn"); - /* dbgassert(free_fn,NULLPTRDEREF "free_fn"); */ - - Workers* w = a_malloc(sizeof(Workers)); - if(!w){ - WARN(ALLOCFAILED "Workers_create"); - return NULL; - } - - w->jq = jq; - w->w_alloc_fcn = init_fn; - w->fn = fn; - w->w_free_fcn = free_fn; - w->global_data = thread_global_data; - w->num_workers = num_workers; - - /* Initialize thread mutex */ - int rv = pthread_mutex_init(&w->global_data_mutex,NULL); - if(rv !=0) { - WARN("Mutex initialization failed"); - return NULL; - } - - - /* Create the threads */ - pthread_t* thread = w->worker_threads; - for(us i = 0; i < num_workers; i++) { - TRACE(15,"Creating thread"); - int rv = pthread_create(thread, - NULL, /* Thread attributes */ - threadfcn, /* Function */ - w); /* Data */ - if(rv!=0) { - WARN("Thread creation failed"); - return NULL; - } - thread++; - } - - return w; - -} - -void Workers_free(Workers* w) { - TRACE(15,"Workers_free"); - dbgassert(w,NULLPTRDEREF "w in Workers_free"); - dbgassert(w->jq,NULLPTRDEREF "w->jq in Workers_free"); - - for(us i=0;inum_workers;i++) { - /* Push the special NULL job. This will make the worker - * threads stop their execution. */ - JobQueue_push(w->jq,NULL); - } - - JobQueue_wait_alldone(w->jq); - - /* Join the threads */ - pthread_t* thread = w->worker_threads; - for(us i=0;inum_workers;i++) { - void* retval; - if(pthread_join(*thread,&retval)!=0) { - WARN("Error joining thread!"); - } - if((retval) != NULL) { - WARN("Thread returned with error status"); - } - thread++; - } - - /* Destroy the global data mutex */ - int rv = pthread_mutex_destroy(&w->global_data_mutex); - if(rv != 0){ - WARN("Mutex destroy failed. Do not know what to do."); - } - - /* All threads joined */ - a_free(w); - -} - -static void* threadfcn(void* thread_global_data) { - - TRACE(15,"Started worker thread function"); - /* dbgassert(thread_global_data,NULLPTRDEREF "thread_data in" */ - /* " threadfcn"); */ - Workers* w = (Workers*) thread_global_data; - - JobQueue* jq = w->jq; - worker_alloc_function walloc = w->w_alloc_fcn; - worker_free_function wfree = w->w_free_fcn; - worker_function worker_fn = w->fn; - void* global_data = w->global_data; - - dbgassert(jq,NULLPTRDEREF "jq in threadfcn"); - /* dbgassert(walloc,NULLPTRDEREF "walloc in threadfcn"); */ - /* dbgassert(wfree,NULLPTRDEREF "wfree in threadfcn"); */ - - int rv = pthread_mutex_lock(&w->global_data_mutex); - if(rv !=0) { - WARN("Global data mutex lock failed"); - pthread_exit((void*) 1); - } - void* w_data = NULL; - if(walloc) { - w_data = walloc(global_data); - if(!w_data) { - WARN(ALLOCFAILED); - pthread_exit((void*) 1); - } - } - - rv = pthread_mutex_unlock(&w->global_data_mutex); - if(rv !=0) { - WARN("Global data mutex unlock failed"); - pthread_exit((void*) 1); - } - - void* job = NULL; - TRACE(20,"Worker ready"); - while (true) { - - TRACE(10,"--------------- START CYCLE -------------"); - job = JobQueue_assign(jq); - - /* Kill the thread for the special NULL job */ - if(!job) break; - - /* Run the worker function */ - rv = worker_fn(w_data,job); - if(rv!=0) { - WARN("An error occured during execution of worker function"); - JobQueue_done(jq,job); - break; - } - - JobQueue_done(jq,job); - TRACE(10,"--------------- CYCLE COMPLETE -------------"); - } - - JobQueue_done(jq,job); - - /* Call the cleanup function */ - if(wfree) { - wfree(w_data); - } - TRACE(15,"Exiting thread. Goodbye"); - pthread_exit((void*) NULL); - - /* This return statement is never reached, but added to have a proper return - * type from this function. */ - return NULL; -} - -#endif // LASP_PARALLEL - -////////////////////////////////////////////////////////////////////// diff --git a/src/lasp/c/lasp_worker.h b/src/lasp/c/lasp_worker.h deleted file mode 100644 index 9780e14..0000000 --- a/src/lasp/c/lasp_worker.h +++ /dev/null @@ -1,65 +0,0 @@ -// worker.h -// -// Author: J.A. de Jong - ASCEE -// -// Description: Provides a clean interface to a pool of worker -// threads. This class is used to easily interface with worker threads -// by just providing a simple function to be called by a worker thread -// on the push of a new job. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef LASP_WORKER_H -#define LASP_WORKER_H -#include "lasp_config.h" -#include "lasp_types.h" - -#ifdef LASP_PARALLEL -typedef struct Workers_s Workers; -typedef struct JobQueue_s JobQueue; - -typedef void* (*worker_alloc_function)(void* global_data); -typedef int (*worker_function)(void* worker_data,void* job); -typedef void (*worker_free_function)(void* worker_data); - -/** - * Create a pool of worker threads, that pull jobs from the job queue - * and perform the action. - * - * @param num_workers Number of worker threads to create - * - * @param jq JobQueue. JobQueue where jobs for the workers are - * pushed. Should stay valid as long as the Workers are alive. - * - * @param worker_alloc_function Function pointer to the function that - * will be called right after the thread has been created. The worker - * alloc function will get a pointer to the thread global data. This - * data will be given to each thread during initialization. Using a - * mutex to avoid race conditions on this global data. - - * @param fn Worker function that performs the action on the - * data. Will be called every time a job is available from the - * JobQueue. Should have a return code of 0 on success. - * - * @param worker_free_function Cleanup function that is called on - * exit. - * - * @return Pointer to Workers handle. NULL on error. - */ -Workers* Workers_create(const us num_workers, - JobQueue* jq, - worker_alloc_function init_fn, - worker_function fn, - worker_free_function free_fn, - void* thread_global_data); - -/** - * Free the pool of workers. - * - * @param w - */ -void Workers_free(Workers* w); - -#endif // LASP_PARALLEL - -#endif // LASP_WORKER_H -//////////////////////////////////////////////////////////////////////