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

Removed old C-code. This is not the way forward anymore

Browse Source
This commit is contained in:
Anne de Jong 2022-09-03 16:10:12 +02:00
parent 34239bfabf
commit 10749137ec
38 changed files with 0 additions and 5255 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
src/lasp/c/CMakeLists.txt
View File

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

128
src/lasp/c/dec_filter_4.c
View File

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

223
src/lasp/c/lasp_alg.c
View File

@ -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;j<A->n_cols;j++){
for(i=0;i<A->n_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 */
//////////////////////////////////////////////////////////////////////

192
src/lasp/c/lasp_alg.h
View File

@ -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;col<y->n_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;col<y->n_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;col<y->n_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;col<y->n_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;col<result->n_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
//////////////////////////////////////////////////////////////////////

33
src/lasp/c/lasp_alloc.h
View File

@ -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 <malloc.h>
#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
//////////////////////////////////////////////////////////////////////

291
src/lasp/c/lasp_aps.c
View File

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

85
src/lasp/c/lasp_aps.h
View File

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

29
src/lasp/c/lasp_assert.c
View File

@ -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 <stdlib.h>
#include <stdio.h>
#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

44
src/lasp/c/lasp_assert.h
View File

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

20
src/lasp/c/lasp_bem.c
View File

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

221
src/lasp/c/lasp_decimation.c
View File

@ -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;channelno<nchannels;channelno++) {
dec->fbs[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;col<downsampled->n_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;chan<nchannels;chan++) {
vd samples_channel = dmat_column((dmat*)samples,
chan);
/* Low-pass filter stuff */
dmat filtered_res = Firfilterbank_filter(dec->fbs[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;chan<dec->nchannels;chan++) {
Firfilterbank_free(dec->fbs[chan]);
}
a_free(dec->fbs);
a_free(dec);
feTRACE(15);
}
//////////////////////////////////////////////////////////////////////

54
src/lasp/c/lasp_decimation.h
View File

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

145
src/lasp/c/lasp_dfifo.c
View File

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

73
src/lasp/c/lasp_dfifo.h
View File

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

70
src/lasp/c/lasp_eq.c
View File

@ -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;filter<eq->nfilters;filter++) {
d ampl = *getvdval(&(eq->ampl_values), filter);
/// TODO: Replace this code with something more fast from BLAS.
for(us sample=0;sample<filtered.n_rows;sample++) {
d* res = getvdval(&result, sample);
*res = *res + *getdmatval(&filtered, sample, filter) * ampl;
}
}
dmat_free(&filtered);
feTRACE(15);
return result;
}
void Eq_setLevels(Eq* eq,const vd* levels) {
fsTRACE(15);
assertvalidptr(eq);
assert_vx(levels);
dbgassert(levels->n_rows == eq->nfilters, "Invalid levels size");
for(us ch=0;ch<eq->nfilters;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);
}

59
src/lasp/c/lasp_eq.h
View File

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

109
src/lasp/c/lasp_fft.h
View File

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

174
src/lasp/c/lasp_firfilterbank.c
View File

@ -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;col<nfilters;col++) {
vc output_fft_col = cmat_column(&output_fft,col);
vc filter_col = cmat_column(&fb->filters,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;
}
//////////////////////////////////////////////////////////////////////

59
src/lasp/c/lasp_firfilterbank.h
View File

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

54
src/lasp/c/lasp_mat.c
View File

@ -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 <math.h>
#if LASP_DEBUG == 1
void print_dmat(const dmat* m) {
fsTRACE(50);
size_t row,col;
for(row=0;row<m->n_rows;row++){
indent_trace();
for(col=0;col<m->n_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;row<m->n_rows;row++){
indent_trace();
for(col=0;col<m->n_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
//////////////////////////////////////////////////////////////////////

597
src/lasp/c/lasp_mat.h
View File

@ -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;col<mat->n_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;col<mat->n_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;col<ncols;col++) {
d* to_d = getdmatval(to,startrow_to,col);
d* from_d = getdmatval(from,startrow_from,col);
d_copy(to_d,from_d,nrows,1,1);
}
}
/**
* 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 dmat dmat_submat(const dmat* parent,
const us startrow,
const us startcol,
const us n_rows,
const us n_cols) {
dbgassert(parent,NULLPTRDEREF);
dbgassert(n_rows+startrow <= parent->n_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;col<to->n_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;col<to->n_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;col<a->n_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;col<x->n_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
//////////////////////////////////////////////////////////////////////

152
src/lasp/c/lasp_math_raw.c
View File

@ -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 <cblas.h>
#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<size;i++) {
res[i] = arr1[i]*arr2[i];
}
#endif
}
void c_hadamard(c res[],
const c arr1[],
const c arr2[],
const us size) {
fsTRACE(15);
uVARTRACE(15,size);
dbgassert(arr1 && arr2 && res,NULLPTRDEREF);
#if LASP_USE_BLAS == 1
#if LASP_DEBUG
if(arr1 == arr2) {
DBGWARN("c_elem_prod_c: 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 /* LASP_DEBUG */
#if LASP_DOUBLE_PRECISION == 1
#define elem_prod_fun cblas_zgbmv
#else
#define elem_prod_fun cblas_cgbmv
#endif
c alpha = 1.0;
c beta = 0.0;
TRACE(15,"Calling " annestr(elem_prod_fun));
uVARTRACE(15,size);
elem_prod_fun(CblasColMajor,
CblasNoTrans,
(blasint) size, /* M: Number of rows */
(blasint) size, /* B: Number of columns */
0, /* KL: Number of sub-diagonals */
0, /* KU: Number of super-diagonals */
(d*) &alpha, /* Multiplication factor */
(d*) arr2, /* A */
1, /* LDA */
(d*) arr1, /* x */
1, /* incX = 1 */
(d*) &beta,
(d*) res, /* The Y matrix to write to */
1); /* incY (increment in res) */
#undef elem_prod_fun
#else /* No blas routines, routine is very simple, but here we
* go! */
DBGWARN("Performing slower non-blas vector-vector multiplication");
for(us i=0;i<size;i++) {
res[i] = arr1[i]*arr2[i];
}
#endif
feTRACE(15);
}
//////////////////////////////////////////////////////////////////////

480
src/lasp/c/lasp_math_raw.h
View File

@ -1,480 +0,0 @@
// lasp_math_raw.h
//
// Author: J.A. de Jong - ASCEE
//
// Description: Raw math routines working on raw arrays of floats and
// complex numbers.
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_MATH_RAW_H
#define LASP_MATH_RAW_H
#ifdef __cplusplus
extern "C" {
#endif
#include "lasp_config.h"
#include "lasp_assert.h"
#include "lasp_tracer.h"
#include "lasp_types.h"
#include <float.h>
#include <math.h>
#if LASP_USE_BLAS == 1
#include <cblas.h>
#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;i<size;i++) {
to[i]=val;
}
}
/**
* Set all elements in an array equal to val
*
* @param to
* @param val
* @param size
*/
static inline void c_set(c to[],c val,us size) {
for(us i=0;i<size;i++) {
to[i]=val;
}
}
/**
* Return the maximum of two doubles
*
* @param a value 1
* @param b value 2
*
* @returns the maximum of value 1 and 2
*/
static inline d d_max(const d a,const d b) {
return a>b?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<size;i++){
result+=a[i]*b[i];
}
return result;
}
/**
* Return the dot product of two complex-valued arrays. Wraps BLAS
* when LASP_USE_BLAS == 1.
*
* @param a complex-valued array
* @param b complex-valued array
* @param size the size of the arrays. *Should be equal-sized!*
*
* @return the dot product
*/
static inline c cc_dot(const c a[],const c b[],us size){
#if LASP_USE_BLAS == 1
WARN("CBlas zdotu not yet tested");
#if LASP_DOUBLE_PRECISION == 1
// Openblas and MKL blas provide some subroutines with a different name, we
// correct for this fact here
#ifdef __MKL_CBLAS_H__
d res;
cblas_zdotu_sub(size,(d*) a,1,(d*) b,1, &res);
return res;
#else
return cblas_zdotu(size,(d*) a,1,(d*) b,1);
#endif
#else
#ifdef __MKL_CBLAS_H__
d res;
cblas_cdotu_sub(size,(d*) a,1,(d*) b,1, &res);
return res;
#endif
return cblas_cdotu(size,(d*) a,1,(d*) b,1);
#endif // LASP_DOUBLE_PRECISION
#else
c result = 0;
us i;
for(i=0;i<size;i++){
result+=a[i]*b[i];
}
return result;
#endif
}
/**
* Compute the dot product of two real arrays.
*
* @param a First array.
* @param b Second array.
* @param size Size of the arrays.
* @return The result.
*/
static inline d d_dot(const d a[],const d b[],const us size){
#if LASP_USE_BLAS == 1
#if LASP_DOUBLE_PRECISION
return cblas_ddot(size,a,1,b,1);
#else // Single precision function
return cblas_sdot(size,a,1,b,1);
#endif
#else // No BLAS, do it manually
d result = 0;
us i;
for(i=0;i<size;i++){
result+=a[i]*b[i];
}
return result;
#endif
}
/**
* Copy array of floats.
*
* @param to : Array to write to
* @param from : Array to read from
* @param size : Size of arrays
* @param to_inc : Mostly equal to 1, the stride of the array to copy to
* @param from_inc : Mostly equal to 1, the stride of the array to copy from
*/
static inline void d_copy(d to[],
const d from[],
const us size,
const us to_inc,
const us from_inc){
#if LASP_USE_BLAS == 1
cblas_dcopy(size,from,from_inc,to,to_inc);
#else
us i;
for(i=0;i<size;i++)
to[i*to_inc] = from[i*from_inc];
#endif
}
/**
* Copy array of floats to array of complex floats. Imaginary part set
* to zero.
*
* @param to : Array to write to
* @param from : Array to read from
* @param size : Size of arrays
*/
static inline void cd_copy(c to[],const d from[],const us size) {
us i;
for(i=0;i<size;i++) {
to[i] = (c) (from[i]);
dbgassert(cimag(to[i]) == 0,"Imaginary part not equal to zero");
}
}
/**
* Copy float vector to complex vector. Imaginary part set
* to zero.
*
* @param to : Vector to write to
* @param from : Vector to read from
*/
static inline void c_copy(c to[],const c from[],const us size){
#if LASP_USE_BLAS == 1
#if LASP_DOUBLE_PRECISION == 1
cblas_zcopy(size,(d*) from,1,(d*) to,1);
#else
cblas_ccopy(size,(d*) from,1,(d*) to,1);
#endif
#else
us i;
for(i=0;i<size;i++)
to[i] = from[i];
#endif
}
/**
* Multiply y with fac, and add result to x
*
* @param x[in,out] Array to add to
* @param y[in] Array to add to x
* @param[in] fac Factor with which to multiply y
* @param[in] size Size of the arrays
*/
static inline void d_add_to(d x[],const d y[],
const d fac,const us size){
#if LASP_USE_BLAS == 1
#if LASP_DOUBLE_PRECISION == 1
cblas_daxpy(size,fac,y,1,x,1);
#else
cblas_saxpy(size,fac,y,1,x,1);
#endif
#else
us i;
for(i=0;i<size;i++)
x[i]+=fac*y[i];
#endif
}
/**
* x = x + fac*y
*
* @param[in,out] x Array to add to
* @param[in] y Array to add to x
* @param[in] fac Factor with which to multiply y
* @param[in] size Size of the arrays
*/
static inline void c_add_to(c x[],const c y[],
const c fac,const us size){
fsTRACE(15);
#if LASP_USE_BLAS == 1
#if LASP_DOUBLE_PRECISION == 1
cblas_zaxpy(size,(d*) &fac,(d*) y,1,(d*) x,1);
#else
cblas_caxpy(size,(d*) &fac,(d*) y,1,(d*) x,1);
#endif
#else
us i;
for(i=0;i<size;i++)
x[i]+=fac*y[i];
#endif
feTRACE(15);
}
/**
* Scale an array of doubles
*
* @param a array
* @param scale_fac scale factor
* @param size size of the array
*/
static inline void d_scale(d a[],const d scale_fac,us size){
#if LASP_USE_BLAS == 1
#if LASP_DOUBLE_PRECISION == 1
cblas_dscal(size,scale_fac,a,1);
#else
cblas_sscal(size,scale_fac,a,1);
#endif
#else
us i;
for(i=0;i<size;i++)
a[i]*=scale_fac;
#endif
}
/**
* Scale an array of complex floats
*
* @param a array
* @param scale_fac scale factor
* @param size size of the array
*/
static inline void c_scale(c a[],const c scale_fac,us size){
#if LASP_USE_BLAS == 1
// Complex argument should be given in as array of two double
// values. The first the real part, the second the imaginary
// part. Fortunately the (c) type stores the two values in this
// order. To be portable and absolutely sure anything goes well,
// we convert it explicitly here.
#if LASP_DOUBLE_PRECISION == 1
cblas_zscal(size,(d*) &scale_fac,(d*) a,1);
#else
cblas_cscal(size,(d*) &scale_fac,(d*) a,1);
#endif
#else
us i;
for(i=0;i<size;i++)
a[i]*=scale_fac;
#endif
}
/**
* Compute the maximum value of an array
*
* @param a array
* @param size size of the array
* @return maximum
*/
static inline d darray_max(const d a[],us size){
us i;
d max = a[0];
for(i=1;i<size;i++){
if(a[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<size;i++){
if(a[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<size;i++){
norm+=a[i]*a[i];
}
norm = d_sqrt(norm);
return norm;
#endif
}
/**
* Compute the \f$ L_2 \f$ norm of an array of complex floats
*
* @param a Array
* @param size Size of array
*/
static inline d c_norm(const c a[],us size){
#if LASP_USE_BLAS == 1
return cblas_dznrm2(size,(d*) a,1);
#else
d norm = 0;
us i;
for(i=0;i<size;i++){
d absa = c_abs(a[i]);
norm+=absa*absa;
}
norm = d_sqrt(norm);
return norm;
#endif
}
/**
* Computes the element-wise vector product, or Hadamard product of
* arr1 and arr2
*
* @param res Where the result will be stored
* @param arr1 Array 1
* @param vec2 Array 2
* @param size: Size of the arrays
*/
void d_elem_prod_d(d res[],
const d arr1[],
const d arr2[],
const us size);
/**
* Computes the element-wise vector product, or Hadamard product of
* arr1 and arr2 for complex floats
*
* @param res Where the result will be stored
* @param arr1 Array 1
* @param vec2 Array 2
* @param size: Size of the arrays
*/
void c_hadamard(c res[],
const c arr1[],
const c arr2[],
const us size);
/**
* Compute the complex conjugate of a complex vector and store the
* result.
*
* @param res Result vector
* @param in Input vector
* @param size Size of the vector
*/
static inline void carray_conj(c res[],const c in[],const us size) {
// First set the result vector to zero
fsTRACE(15);
c_set(res,0,size);
#if LASP_USE_BLAS == 1
#if LASP_DOUBLE_PRECISION == 1
// Cast as a float, scale all odd elements with minus one to find
// the complex conjugate.
cblas_daxpy(size ,1.0,(d*) in,2,(d*) res,2);
cblas_daxpy(size,-1.0,&((d*) in)[1],2,&((d*) res)[1],2);
#else
cblas_faxpy(size ,1,(d*) in,2,(d*) res,2);
cblas_faxpy(size,-1,&((d*) in)[1],2,&((d*) res)[1],2);
#endif // LASP_DOUBLE_PRECISION
#else
for(us i=0;i<size;i++) {
res[i] = c_conj(in[i]);
}
#endif // LASP_USE_BLAS
feTRACE(15);
}
/**
* In place complex conjugation
*
* @param res Result vector
* @param size Size of the vector
*/
static inline void c_conj_inplace(c res[],us size) {
#if LASP_USE_BLAS
#if LASP_DOUBLE_PRECISION == 1
// Cast as a float, scale all odd elements with minus one to find
// the complex conjugate.
cblas_dscal(size,-1,&((d*) res)[1],2);
#else
cblas_sscal(size,-1,&((d*) res)[1],2);
#endif // LASP_DOUBLE_PRECISION
#else
for(us i=0;i<size;i++) {
res[i] = c_conj(res[i]);
}
#endif // LASP_USE_BLAS
}
#ifdef __cplusplus
}
#endif
#endif // LASP_MATH_RAW_H
//////////////////////////////////////////////////////////////////////

367
src/lasp/c/lasp_mq.c
View File

@ -1,367 +0,0 @@
// lasp_mq.c
//
// Author: J.A. de Jong -ASCEE
//
// Description: Message queue implementation using a linked
// list. Using mutexes and condition variables to sync access between
// different threads.
//////////////////////////////////////////////////////////////////////
#define TRACERPLUS (-6)
#include "lasp_types.h"
#include "lasp_tracer.h"
#include "lasp_assert.h"
#include "lasp_alloc.h"
#include "lasp_mq.h"
#include <pthread.h>
/* #ifdef linux */
#define _GNU_SOURCE
#include <sys/types.h>
#include <unistd.h>
#ifndef MS_WIN64
#include <sys/syscall.h>
#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;i<jq->max_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;i<jq->max_jobs;i++){
if(j->ready && !j->running)
return j;
j++;
}
return NULL;
}
void print_job_queue(JobQueue* jq) {
fsTRACE(15);
for(us i=0;i<jq->max_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;jindex<max_jobs;jindex++) {
j->job_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;i<max_jobs;i++) {
if(j->ready == 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;i<jq->max_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
//////////////////////////////////////////////////////////////////////

78
src/lasp/c/lasp_mq.h
View File

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

18
src/lasp/c/lasp_nprocs.c
View File

@ -1,18 +0,0 @@
#include "lasp_nprocs.h"
#ifdef MS_WIN64
#include <windows.h>
#else
// Used for obtaining the number of processors
#include <sys/sysinfo.h>
#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
}

25
src/lasp/c/lasp_nprocs.h
View File

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

181
src/lasp/c/lasp_ps.c
View File

@ -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;i<nchannels;i++) {
vd column = dmat_column((dmat*) timedata,i);
vd column_work = dmat_column(&timedata_work,i);
vd_elem_prod(&column_work,&column,&ps->window);
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<nchannels;i++) {
for(us j=0;j<nchannels;j++) {
/* The indices here are important. This is also how it
* is documented */
vc res = cmat_column(result,i+j*nchannels);
vc i_vec = cmat_column(&fft_work,i);
vc j_vec = cmat_column(&fft_work,j);
/* Compute the conjugate of spectra j */
cmat_conj(&j_vec_conj,&j_vec);
check_overflow_xmat(fft_work);
/* Compute the element-wise product of the two vectors and
* store the result as the result */
vc_hadamard(&res,&i_vec,&j_vec_conj);
vc_free(&i_vec);
vc_free(&j_vec);
vc_free(&res);
}
}
check_overflow_xmat(*result);
check_overflow_xmat(*timedata);
cmat_free(&fft_work);
vc_free(&j_vec_conj);
feTRACE(15);
}
//////////////////////////////////////////////////////////////////////

61
src/lasp/c/lasp_ps.h
View File

@ -1,61 +0,0 @@
// lasp_ps.h
//
// Author: J.A. de Jong - ASCEE
//
// Description: Single sided power and cross-power spectra computation
// routines.
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_PS_H
#define LASP_PS_H
#include "lasp_window.h"
typedef struct PowerSpectra_s PowerSpectra;
/**
* Allocate a PowerSpectra computer
*
* @param nfft fft length
* @param nchannels Number of channels
* @param wt Windowtype, as defined in window.h
*
* @return PowerSpectra handle, NULL on error
*/
PowerSpectra* PowerSpectra_alloc(const us nfft,
const WindowType wt);
/**
* Compute power spectra.
*
* @param ps[in] PowerSpectra handle
* @param[in] timedata Time data. Should have size nfft*nchannels
*
* @param[out] result Here, the result will be stored. Should have
* size (nfft/2+1)*nchannels^2, such that the Cij at frequency index f
* can be obtained as result[f,i+j*nchannels]
*
*/
void PowerSpectra_compute(const PowerSpectra* ps,
const dmat* timedata,
cmat *result);
/**
* Return nfft
*
* @param ps[in] PowerSpectra handle
*
* @return nfft
*/
us PowerSpectra_getnfft(const PowerSpectra* ps);
/**
* Free PowerSpectra
*
* @param[in] ps PowerSpectra handle
*/
void PowerSpectra_free(PowerSpectra* ps);
#endif // LASP_PS_H
//////////////////////////////////////////////////////////////////////

108
src/lasp/c/lasp_pyarray.h
View File

@ -1,108 +0,0 @@
// ascee_python.h
//
// Author: J.A. de Jong - ASCEE
//
// Description:
// Routine to generate a Numpy array from an arbitrary buffer. Careful, this
// code should both be C and C++ compatible!!
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_PYARRAY_H
#define LASP_PYARRAY_H
#include "lasp_types.h"
#include <assert.h>
#include <numpy/ndarrayobject.h>
#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 <malloc.h>
/**
* 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

59
src/lasp/c/lasp_python.h
View File

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

32
src/lasp/c/lasp_signals.h
View File

@ -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;i<signal->n_rows;i++) {
res+= d_pow(*getvdval(signal,i),2);
}
res /= signal->n_rows;
return res;
}
#endif // LASP_SIGNALS_H
//////////////////////////////////////////////////////////////////////

293
src/lasp/c/lasp_slm.c
View File

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

96
src/lasp/c/lasp_slm.h
View File

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

264
src/lasp/c/lasp_sosfilterbank.c
View File

@ -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 a<b?a:b;
}
static inline us max(us a, us b){
return a>b?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;index<nsections*6;index++){
// Copy contents to position in sos matrix
*getdmatval(sos,filter_no,index) = *getvdval(&filter_coefs,index);
}
feTRACE(15);
}
void Sosfilterbank_free(Sosfilterbank* fb) {
fsTRACE(15);
assertvalidptr(fb);
dmat_free(&(fb->sos));
dmat_free(&(fb->state));
#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;section<nsections;section++) {
d w1 = *getdmatval(&(job->state),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;sample<nsamples;sample++){
d w0 = *y - a1*w1 - a2*w2;
d yn = b0*w0 + b1*w1 + b2*w2;
w2 = w1;
w1 = w0;
*y++ = yn;
}
*getdmatval(&(job->state),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;filter<filterbank_size;filter++) {
d_copy(getdmatval(&ys,0,filter),getvdval(xs,0),nsamples,1,1);
}
Job jobs[MAX_SOS_FILTER_BANK_SIZE];
/// Implementation is based on Proakis & Manolakis - Digital Signal
/// Processing, Fourth Edition, p. 550
Job job_template = {0, nsections, nsamples, sos, state, ys};
for(us filter=0;filter<filterbank_size;filter++) {
/// Obtain state information for current section, and all filters
jobs[filter] = job_template;
jobs[filter].filter_no = filter;
#ifdef LASP_PARALLEL
if(fb->workers) {
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;
}

81
src/lasp/c/lasp_sosfilterbank.h
View File

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

206
src/lasp/c/lasp_worker.c
View File

@ -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 <pthread.h>
#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;i<w->num_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;i<w->num_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
//////////////////////////////////////////////////////////////////////

65
src/lasp/c/lasp_worker.h
View File

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