commit 52e5a31bdbf6eb1c2bafbbc071d9feb5d7fe3033 Author: J.A. de Jong Date: Mon Jan 29 16:14:50 2018 +0100 Initial commit. Lots of stuff diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b5e9a21 --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +*.a +CMakeFiles +CMakeCache.txt +cmake_install.cmake +Beamforming.egg-info +beamforming/*.cpp +Makefile +build +*.html +__pycache__ +cython_debug +beamforming/config.pxi +beamforming/fft.c +*.so diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..4546138 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,134 @@ +cmake_minimum_required (VERSION 3.0) +project(beamforming) + +# Whether we want to use blas yes or no +set(ASCEE_USE_BLAS 1) +add_definitions(-DASCEE_USE_BLAS=1) + +set(ASCEE_FLOAT double) +# set(ASCEE_FLOAT float) + +add_definitions(-DASCEE_PARALLEL) +add_definitions(-DASCEE_MAX_NUM_THREADS=8) + +add_definitions(-DASCEE_MAX_NUM_CHANNELS=80) + +# ####################################### End of user-adjustable variables section + +add_definitions(-D_REENTRANT) + + +if(ASCEE_FLOAT STREQUAL "double") + add_definitions(-DASCEE_FLOAT=64) + add_definitions(-DASCEE_DOUBLE_PRECISION) +else() + add_definitions(-DASCEE_FLOAT=32) +endif(ASCEE_FLOAT STREQUAL "double") + + +if(NOT DEFINED ASCEE_DEBUG) + message(SEND_ERROR "ASCEE_DEBUG flag not defined. Please set -DASCEE_DEBUG=TRUE or -DASCEE_DEBUG=FALSE") +endif(NOT DEFINED ASCEE_DEBUG) + +# ##################### END Cmake variables converted to a macro +set(Python_ADDITIONAL_VERSIONS "3") +# #################### Setting definitions and debug-specific compilation flags + +if(ASCEE_DEBUG) + set(TRACERNAME ASCEETracer) + set(CMAKE_BUILD_TYPE Debug) + message("Building debug code") + set(CMAKE_BUILD_TYPE Debug) + add_definitions(-DASCEE_DEBUG=1) + add_definitions(-DTRACERNAME=${TRACERNAME}) + add_definitions(-DDEBUG) + add_definitions(-DTRACER=1) + + set(CYTHON_VARIABLES "#cython: boundscheck=True") + # This will produce html files + set(CYTHON_ANNOTATE ON) + # Add the __FILENAME__ macro + # set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -D__FILENAME__='\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"'") +else() + message("Building release code") + set(CMAKE_BUILD_TYPE Release) + set(CYTHON_ANNOTATE OFF) + set(CYTHON_NO_DOCSTRINGS ON) + set(CYTHON_VARIABLES "#cython: boundscheck=False,wraparound=False,embedsignature=False,emit_code_comments=False") + # Strip unnecessary symbols + # set(CMAKE_SHARED_LINKER_FLAGS "-Wl,--gc-sections") + # set(CMAKE_MODULE_LINKER_FLAGS "-Wl,--gc-sections") + + add_definitions(-DTRACER=0 -DNDEBUG) +endif(ASCEE_DEBUG) + +# The last argument here takes care of calling SIGABRT when an integer overflow +# occures. +############################## General compilation flags (independent of debug mode, windows or linux) + +# General make flags +set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fPIC -std=c11 -Wall -Wextra \ + -Wimplicit-function-declaration -ftrapv -Wno-type-limits") + +# Debug make flags +set(CMAKE_C_FLAGS_DEBUG "-g" ) + +set(CMAKE_C_FLAGS_RELEASE "-O2 -mfpmath=sse -march=x86-64 -mtune=native \ + -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions") + +# set(CMAKE_C_FLAGS_RELEASE "-O2 -march=native -mtune=native -fomit-frame-pointer") + + +if(ASCEE_USE_BLAS) + add_definitions(-DASCEE_USE_BLAS=1) +else() + add_definitions(-DASCEE_USE_BLAS=0) +endif(ASCEE_USE_BLAS) + +if(CMAKE_SYSTEM_NAME STREQUAL "Windows") + set(win32 true) +else() + set(win32 false) +endif(CMAKE_SYSTEM_NAME STREQUAL "Windows") + +# If numpy cannot be found in the standard include path of the Python +if(DEFINED NUMPY_INCLUDE) + include_directories(${NUMPY_INCLUDE}) +endif(DEFINED NUMPY_INCLUDE) + + + + +# ############################# End compilation flags + +add_subdirectory(fftpack) +include_directories( + fftpack + beamforming/c + ) +add_subdirectory(beamforming) +add_subdirectory(test) + +find_program(PYTHON "python") + +if (PYTHON) + # set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in") + set(SETUP_PY "${CMAKE_CURRENT_BINARY_DIR}/setup.py") + set(DEPS "${CMAKE_CURRENT_SOURCE_DIR}/beamforming/*.py") + set(OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp") + + # configure_file(${SETUP_PY_IN} ${SETUP_PY}) + + add_custom_command(OUTPUT ${OUTPUT} + COMMAND ${PYTHON} ${SETUP_PY} build + COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT} + DEPENDS ${DEPS}) + + add_custom_target(target ALL DEPENDS ${OUTPUT}) + + install(CODE "execute_process(COMMAND ${PYTHON} ${SETUP_PY} install)") +endif() + + +############################## End compiler settings + diff --git a/beamforming/CMakeLists.txt b/beamforming/CMakeLists.txt new file mode 100644 index 0000000..2311f17 --- /dev/null +++ b/beamforming/CMakeLists.txt @@ -0,0 +1,19 @@ +add_subdirectory(c) + +configure_file(config.pxi.in config.pxi) + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/cmake") +find_package( PythonLibs REQUIRED ) +include(UseCython) + + +include_directories( + . + c + ) + + +set_source_files_properties(fft.c PROPERTIES COMPILE_FLAGS + "-Wno-cpp -Wno-sign-compare -Wno-incompatible-pointer-types -Wno-strict-aliasing -Wno-implicit-fallthrough") +cython_add_module(fft fft.pyx) +target_link_libraries(fft beamforming_lib ) diff --git a/beamforming/__init__.py b/beamforming/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/beamforming/c/CMakeLists.txt b/beamforming/c/CMakeLists.txt new file mode 100644 index 0000000..d33e7c1 --- /dev/null +++ b/beamforming/c/CMakeLists.txt @@ -0,0 +1,18 @@ +if(!ASCEE_DEBUG) +# SET_SOURCE_FILES_PROPERTIES(si_lpn.c PROPERTIES COMPILE_FLAGS -O3) +# SET_SOURCE_FILES_PROPERTIES(si_math.c PROPERTIES COMPILE_FLAGS -O3) +endif(!ASCEE_DEBUG) + +add_library(beamforming_lib + fft.c + ascee_math.c + ascee_assert.c + tracer.c + window.c + # ps.c + mq.c + worker.c + ) + + +target_link_libraries(beamforming_lib fftpack openblas) diff --git a/beamforming/c/ascee_alloc.h b/beamforming/c/ascee_alloc.h new file mode 100644 index 0000000..4ec31ad --- /dev/null +++ b/beamforming/c/ascee_alloc.h @@ -0,0 +1,22 @@ +// ascee_alloc.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// memory allocation functions. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef ASCEE_ALLOC_H +#define ASCEE_ALLOC_H +#include + +/** + * Reserved words for memory allocation. Can be changed to something + * else when required. For example for debugging purposes. + */ +#define a_malloc malloc +#define a_free free +#define a_realloc realloc + +#endif // ASCEE_ALLOC_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/ascee_assert.c b/beamforming/c/ascee_assert.c new file mode 100644 index 0000000..6eacb11 --- /dev/null +++ b/beamforming/c/ascee_assert.c @@ -0,0 +1,28 @@ +// ascee_assert.c +// +// Author: J.A. de Jong -ASCEE +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#include "ascee_assert.h" + +#ifdef ASCEE_DEBUG +#include +#include +#include "types.h" +#define MAX_MSG 200 + +void DBG_AssertFailedExtImplementation(const char* filename, + us linenr, + const char* extendedinfo) +{ + char scratchpad[MAX_MSG]; + + sprintf(scratchpad,"ASSERT: file %s line %lu: (%s)\n", + filename, linenr, extendedinfo); + printf("%s\n", scratchpad); + abort(); +} + +#endif diff --git a/beamforming/c/ascee_assert.h b/beamforming/c/ascee_assert.h new file mode 100644 index 0000000..46e134c --- /dev/null +++ b/beamforming/c/ascee_assert.h @@ -0,0 +1,40 @@ +// ascee_assert.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// Basic tools for debugging using assert statements including text. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef ASCEE_ASSERT_H +#define ASCEE_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 ASCEE_DEBUG +#include "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 ); \ + } + +#else // ASCEE_DEBUG not defined + +#define dbgassert(assertion, assert_string) + +#endif // ASCEE_DEBUG + +#endif // ASCEE_ASSERT_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/ascee_math.c b/beamforming/c/ascee_math.c new file mode 100644 index 0000000..fa53f2d --- /dev/null +++ b/beamforming/c/ascee_math.c @@ -0,0 +1,453 @@ +// si_math.c +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#define TRACERPLUS (-10) + +#include "ascee_assert.h" +#include "ascee_math.h" +#include "tracer.h" + +#if ASCEE_USE_BLAS +#include +#endif + +#include + +#ifdef ASCEE_DEBUG +void print_cmat(const cmat* m) { + size_t row,col; + for(row=0;rown_rows;row++){ + for(col=0;coln_cols;col++){ + c val = m->data[row+m->n_rows*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"); + + } +} +void print_vc(const vc* m) { + TRACE(20,"print_vc"); + size_t row; + + for(row=0;rowsize;row++){ + + d rval = creal(m->data[row]); + d ival = cimag(m->data[row]); + + printf("%c%2.2e%c%2.2ei ",rval< 0 ?'-': ' ', d_abs(rval),ival<0 ? '-' : '+',d_abs(ival) ) ; + printf("\n"); + + } +} +void print_vd(const vd* m) { + TRACE(20,"print_vd"); + size_t row; + iVARTRACE(20,m->size); + for(row=0;rowsize;row++){ + + d rval = m->data[row]; + + printf("%c%2.2e ",rval< 0 ? '\r': ' ',rval); + printf("\n"); + } +} +void print_dmat(const dmat* m) { + size_t row,col; + for(row=0;rown_rows;row++){ + for(col=0;coln_cols;col++){ + d val = m->data[row+m->n_rows*col]; + printf("%c%2.2e ", val<0?'-':' ' ,d_abs(val)); + + } + printf("\n"); + + } +} +#endif + +void d_elem_prod_d(d res[], + const d arr1[], + const d arr2[], + const us size) { + + #if ASCEE_USE_BLAS + + #if ASCEE_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 ASCEE_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;in_rows == b->size); + assert(A->n_cols == x->size); + + #if ASCEE_USE_BLAS == 1 + + /* 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; + size_t n_rows = A->n_rows; + + vc_set(b,0.0); + + iVARTRACE(20,A->n_cols); + iVARTRACE(20,A->n_rows); + + for(j=0;jn_cols;j++){ + for(i=0;in_rows;i++) { + + c* Aij = &A->data[i+j*n_rows]; + b->data[i] += *Aij * x->data[j]; + + } + + } + + + #endif +} + +void kronecker_product(const cmat* a,const cmat* b,cmat* result){ + + assert(result->n_rows == a->n_rows*b->n_rows); + assert(result->n_cols == a->n_cols*b->n_cols); + + c a_rs; + c b_vw; + + int r_col; + int r_row; + + for(size_t r=0; r< a->n_rows;r++){ + + for(size_t s=0; s n_cols;s++) { + + for(size_t v=0;v < b->n_rows; v++) { + + for(size_t w=0;w < b->n_cols;w++) { + + a_rs = *getcmatval(a,r,s); + b_vw = *getcmatval(b,v,w); + + r_row = b->n_rows*r+v; + r_col = b->n_cols*s+w; + + result->data[r_row + r_col * result->n_rows] = a_rs * b_vw; + + } + } + } + } +} /* void kronecker_product */ + +/* #include */ +/* 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 ASCEE_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){ */ + +/* POOL_INIT(lsq_solve_pool); */ +/* 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; */ +/* } */ + +/* Pool_free(&lsq_solve_pool,A_data); */ +/* Pool_free(&lsq_solve_pool,work_data); */ +/* POOL_EXIT(lsq_solve_pool,15); */ +/* 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); */ +/* } */ + +////////////////////////////////////////////////////////////////////// + diff --git a/beamforming/c/ascee_math.h b/beamforming/c/ascee_math.h new file mode 100644 index 0000000..0ad39e1 --- /dev/null +++ b/beamforming/c/ascee_math.h @@ -0,0 +1,752 @@ +// ascee_math.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef ASCEE_MATH_H +#define ASCEE_MATH_H +#include +#include "types.h" +#include "tracer.h" +#include "ascee_assert.h" +#include "ascee_alloc.h" + +#if ASCEE_USE_BLAS == 1 +#include +#endif + +#ifdef ASCEE_DOUBLE_PRECISION +#define c_real creal +#define c_imag cimag +#define d_abs fabs +#define c_abs cabs +#define c_conj conj +#define d_atan2 atan2 +#define d_acos acos +#define d_sqrt sqrt +#define c_exp cexp +#define d_sin sin +#define d_cos cos +#define d_pow pow + +#else // ASCEE_DOUBLE_PRECISION not defined +#define c_conj conjf +#define c_real crealf +#define c_imag cimagf +#define d_abs fabsf +#define c_abs cabsf +#define d_atan2 atan2f +#define d_acos acosf +#define d_sqrt sqrtf +#define c_exp cexpf +#define d_sin sinf +#define d_cos cosf +#define d_pow powf + +#endif // ASCEE_DOUBLE_PRECISION + +#ifdef M_PI +static const d number_pi = M_PI; +#else +static const d number_pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679; +#endif + + +/// Code generation for vector of floats and vector of complex floats. +#define vxinit(type) \ + typedef struct { \ + us size; \ + type* data; \ + } v##type; +vxinit(d); +vxinit(c); + +/// Code generation for matrix of floats and matrix of complex floats. +#define xmatinit(type) \ + typedef struct { \ + us n_rows; \ + us n_cols; \ + type* data; \ + } type##mat; +xmatinit(d); +xmatinit(c); + + + +#define setvecval(vec,index,val) \ + dbgassert((((us) index) <= (vec)->size),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)->n_rows+(row)] = val; + +/** + * Return a value from a vector + * + * @param mat The vector + * @param row The row + */ +static inline d* getvdval(const vd* vec,us row){ + dbgassert(row < vec->size,OUTOFBOUNDSVEC); + return &vec->data[row]; +} + +/** + * Return a value from a complex vector + * + * @param mat The vector + * @param row The row + */ +static inline c* getvcval(const vc* vec,us row){ + dbgassert(row < vec->size,OUTOFBOUNDSVEC); + return &vec->data[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){ + assert((row) < mat->n_rows); + assert((col) < mat->n_cols); + return &mat->data[(col)*mat->n_rows+(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(row < mat->n_rows,OUTOFBOUNDSMATR); + dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC); + return &mat->data[col*mat->n_rows+row]; +} + +/** + * Sets all values in a vector to the value + * + * @param b the vector to set + * @param value + */ +static inline void vd_set(vd* vec, d value){ + us i; + for(i=0;isize;i++){ + vec->data[i] = value; + } +} + +/** + * Sets all values in a vector to the value + * + * @param vec the vector to set + * @param value + */ +static inline void vc_set(vc* vec,const c value){ + us i; + for(i=0;isize;i++){ + vec->data[i] = value; + } +} + +/** + * 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){ + us i,size = mat->n_cols*mat->n_rows; + for(i=0;idata[i] = value; + } +} + + +/** + * 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){ + us i,size = mat->n_cols*mat->n_rows; + for(i=0;idata[i] = value; + } +} + +/** + * Return a column pointer of the matrix + * + * @param mtrx The matrix. + * @param column The column number. + * + * @return Pointer to the column. + */ +static inline d* d_column(dmat* mtrx,us column){ + return &mtrx->data[mtrx->n_rows*column]; +} + +/** + * Return a column pointer of the matrix + * + * @param mtrx The matrix. + * @param column The column number. + * + * @return Pointer to the column. + */ +static inline c* c_column(cmat* mtrx,us column){ + return &mtrx->data[mtrx->n_rows*column]; +} + + +/** + * 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 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;isize == b->size,SIZEINEQUAL); + return d_dot(a->data,b->data,a->size); +} + + +/** + * Copy array of floats. + * + * @param to : Array to write to + * @param from : Array to read from + * @param size : Size of arrays + */ +static inline void d_copy(d to[],const d from[],const us size){ + #if ASCEE_USE_BLAS == 1 + cblas_dcopy(size,from,1,to,1); + #else + us i; + for(i=0;isize==from->size,SIZEINEQUAL); + d_copy(to->data,from->data,to->size); +} + +/** + * 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 max) max=a[i]; + } + return max; +} +/** + * Compute the minimum of an array + * + * @param a array + * @param size size of the array + * @return minimum + */ +static inline d d_min(const d a[],us size){ + us i; + d min = a[0]; + for(i=1;i min) min=a[i]; + } + return min; +} + +/** + * Compute the \f$ L_2 \f$ norm of an array of doubles + * + * @param a Array + * @param size Size of array + */ +static inline d d_norm(const d a[],us size){ + #if ASCEE_USE_BLAS == 1 + return cblas_dnrm2(size,a,1); + #else + d norm = 0; + us i; + for(i=0;idata); \ + } +matvec_free(vd); +matvec_free(vc); +matvec_free(dmat); +matvec_free(cmat); + +/** + * Now the following functions exist: vd_free, vc_free, dmat_free and + * cmat_free. + */ + + +/** + * 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, NULL}; + result.data = (d*) a_malloc(n_rows*n_cols*sizeof(d)); + #ifdef ASCEE_DEBUG + dmat_set(&result,NAN); + #endif // ASCEE_DEBUG + assert(result.data); + return result; +} + + +/** + * Allocate data for 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(us n_rows, + us n_cols) { + cmat result = { n_rows, n_cols, NULL}; + result.data = (c*) a_malloc(n_rows*n_cols*sizeof(c)); + #ifdef ASCEE_DEBUG + cmat_set(&result,NAN+I*NAN); + #endif // ASCEE_DEBUG + assert(result.data); + return result; +} + +/** + * Resize an existing dmat or a cmat + */ +#define type_mat_resize(type) \ + static inline void type##mat_resize(type##mat * mat,\ + us nrows,us ncols) { \ + mat->n_rows = nrows; \ + mat->n_cols = ncols; \ + mat->data = realloc(mat->data,(nrows*ncols)*sizeof( type )); \ + } +type_mat_resize(d); +type_mat_resize(c); + +/** + * 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 copy_dmat_rows(dmat* to,const dmat* from, + us startrow_from, + us startrow_to, + us nrows) { + us col,ncols = to->n_cols; + + dbgassert(startrow_from+nrows <= from->n_rows,OUTOFBOUNDSMATR); + dbgassert(startrow_to+nrows <= to->n_rows,OUTOFBOUNDSMATR); + for(col=0;col ASCEE_MAX_NUM_CHANNELS) { + WARN("Too high number of channels! Please increase the " + "ASCEE_MAX_NUM_CHANNELS compilation flag"); + return NULL; + } + + Fft* fft = a_malloc(sizeof(Fft)); + if(fft==NULL) { + WARN("Fft allocation failed"); + return NULL; + } + + fft->nfft = nfft; + fft->nchannels = nchannels; + + fft->fftr = Fftr_alloc(nfft); + if(!fft->fftr) { + WARN(ALLOCFAILED "fftr"); + return NULL; + } + #endif // ASCEE_PARALLEL + + feTRACE(15); + return fft; +} +void Fft_free(Fft* fft) { + fsTRACE(15); + Fftr_free(fft->fftr); + a_free(fft); + feTRACE(15); +} +us Fft_nchannels(const Fft* fft) {return fft->nchannels;} +us Fft_nfft(const Fft* fft) {return fft->nfft;} + +void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result) { + fsTRACE(15); + + dbgassert(timedata->n_rows == fft->nfft,"Invalid size for time data rows." + " Should be equal to nfft"); + dbgassert(timedata->n_cols == fft->nchannels,"Invalid size for time data cols." + " Should be equal to nchannels"); + + for(us col=0;colnchannels;col++) { + + Fftr_fftr(fft->fftr, + getdmatval(timedata,0,col), + getcmatval(result,0,col)); + + } + + feTRACE(15); +} + +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/fft.h b/beamforming/c/fft.h new file mode 100644 index 0000000..96d357c --- /dev/null +++ b/beamforming/c/fft.h @@ -0,0 +1,59 @@ +// fft.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// Interface to the FFT library, multiple channel FFT's +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef FFT_H +#define FFT_H +#include "types.h" +#include "ascee_math.h" + +/** + * Perform forward FFT's on real time data. + * + */ +typedef struct Fft_s Fft; + +/** + * Construct an Fft object + * + * @param nfft Nfft size + * @param nchannels Number of channels + * + * @return Pointer to Fft handle, NULL on error + */ +Fft* Fft_alloc(const us nfft,const us nchannels); + +/** + * Returns the number of channels for this Fft instance + * + * @return nchannels + */ +us Fft_nchannels(const Fft*); + +/** + * Returns the nfft for this Fft instance + * + * @return nfft + */ +us Fft_nfft(const Fft*); + + +/** + * Compute the fft of the data matrix, first axis is assumed to be + * the time axis. + * + * @param[in] timedata Input time data pointer, such that + * data[i*nfft+j) is the i-th channel from data stream j. + * + * @param[out] result: Fft't data, size should be (nfft/2+1)*nchannels + */ +void Fft_fft(const Fft*,const dmat* timedata,cmat* result); + +void Fft_free(Fft* fft); + +#endif // FFT_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/mq.c b/beamforming/c/mq.c new file mode 100644 index 0000000..5de38de --- /dev/null +++ b/beamforming/c/mq.c @@ -0,0 +1,339 @@ +// 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 "types.h" +#include "tracer.h" +#include "ascee_assert.h" +#include "ascee_alloc.h" +#include "mq.h" +#include + +#ifdef __linux__ +#include +#include +#include +#endif + +typedef struct { + void* job_ptr; + bool running; + bool ready; +} Job; + +typedef struct JobQueue_s { + 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. */ + + Job* jobs; /**< Pointer to job vector */ + us max_jobs; /**< Stores the maximum number of + * items */ +} JobQueue; + +static us count_jobs(JobQueue* jq) { + fsTRACE(15); + us njobs = 0; + for(us i=0;imax_jobs;i++){ + if(jq->jobs[i].ready) + njobs++; + } + return njobs; +} +static Job* get_ready_job(JobQueue* jq) { + fsTRACE(15); + Job* j = jq->jobs; + for(us i=0;imax_jobs;i++){ + if(j->ready && !j->running) + return j; + j++; + } + return NULL; +} +void print_job_queue(JobQueue* jq) { + fsTRACE(15); + for(us i=0;imax_jobs;i++) { + printf("Job %zu", i); + if(jq->jobs[i].ready) + printf(" available"); + if(jq->jobs[i].running) + printf(" running"); + + printf(" - ptr %zu\n", (us) jq->jobs[i].job_ptr); + + } + feTRACE(15); +} + + +#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"); \ + } + +JobQueue* JobQueue_alloc(const us max_jobs) { + TRACE(15,"JobQueue_alloc"); + if(max_jobs > ASCEE_MAX_NUM_CHANNELS) { + WARN("Max jobs restricted to ASCEE_MAX_NUM_CHANNELS"); + return NULL; + } + JobQueue* jq = a_malloc(sizeof(JobQueue)); + + + if(!jq) { + WARN("Allocation of JobQueue failed"); + return NULL; + } + jq->max_jobs = max_jobs; + + jq->jobs = a_malloc(max_jobs*sizeof(Job)); + if(!jq->jobs) { + WARN("Allocation of JobQueue jobs failed"); + return NULL; + } + + Job* j = jq->jobs; + for(us jindex=0;jindexjob_ptr = NULL; + j->ready = false; + j->running = false; + j++; + } + + /* 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; + } + + /* 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); + + /* 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."); + } + +} + +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; + + /* 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"); + } + } + + dbgassert(count_jobs(jq) != max_jobs, + "Queue cannot be full!"); + + /* Queue is not full try to find a place, fill it */ + Job* j = jq->jobs; + us i; + for(i=0;iready == false ) { + dbgassert(j->job_ptr==NULL,"Job ptr should be 0"); + dbgassert(j->ready==false,"Job cannot be assigned"); + break; + } + j++; + } + dbgassert(i!=jq->max_jobs,"Should have found a job!"); + + j->job_ptr = job_ptr; + j->ready = true; + + /* 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"); + } + + } + + /* print_job_queue(jq); */ + + UNLOCK_MUTEX; + + return SUCCESS; +} +void* JobQueue_assign(JobQueue* jq) { + + TRACE(15,"JobQueue_assign"); + + LOCK_MUTEX; + + /* Wait until a job is available */ + Job* j; + while ((j=get_ready_job(jq))==NULL) { + + TRACE(15,"JobQueue_assign: no ready job"); + pthread_cond_wait(&jq->cv_plus,&jq->mutex); + + } + + TRACE(16,"JobQueue_assign: found ready job. Assigned to:"); + #ifdef ASCEE_DEBUG + #ifdef __linux__ + + pid_t tid = syscall(SYS_gettid); + iVARTRACE(16,tid); + #endif // __linux__ + + + #endif // ASCEE_DEBUG + + + /* print_job_queue(jq); */ + /* Find a job from the queue, assign it and return it */ + j->running = true; + + 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"); + } + } + + UNLOCK_MUTEX; + + TRACE(15,"End JobQueue_assign"); + + return j->job_ptr; +} +void JobQueue_done(JobQueue* jq,void* job_ptr) { + + TRACE(15,"JobQueue_done"); + dbgassert(jq,NULLPTRDEREF "jq in JobQueue_done"); + + LOCK_MUTEX; + + /* print_job_queue(jq); */ + + /* Find the job from the queue, belonging to the job_ptr */ + Job* j=jq->jobs; + us i; + for(i=0;imax_jobs;i++) { + iVARTRACE(10,i); + if(j->ready && j->running && j->job_ptr == job_ptr) { + TRACE(15,"Found the job that has been done:"); + j->ready = false; + j->job_ptr = NULL; + j->running = false; + break; + } + j++; + } + + /* print_job_queue(jq); */ + + /* Job done, broadcast this */ + rv = pthread_cond_signal(&jq->cv_minus); + if(rv !=0) { + WARN("Condition variable broadcast failed"); + } + + UNLOCK_MUTEX; +} + +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; + +} + + +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/mq.h b/beamforming/c/mq.h new file mode 100644 index 0000000..7feff96 --- /dev/null +++ b/beamforming/c/mq.h @@ -0,0 +1,73 @@ +// mq.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// Multithreaded job queue implementation +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef MQ_H +#define MQ_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 indefinetely until some job is + * available. + * + * @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* jp,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. + * + */ +void JobQueue_wait_alldone(JobQueue*); + +#endif // MQ_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/ps.c b/beamforming/c/ps.c new file mode 100644 index 0000000..8f330d0 --- /dev/null +++ b/beamforming/c/ps.c @@ -0,0 +1,304 @@ +// ps.c +// +// Author: J.A. de Jong -ASCEE +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#define TRACERPLUS 1000 +#include "ps.h" +#include "fft.h" +#include "ascee_alloc.h" +#include "ascee_math.h" +#include "ascee_assert.h" + +typedef struct PowerSpectra_s { + + vd window; + + d win_pow; /**< The power of the window */ + Fft* fft; /**< Fft routines storage */ + cmat fft_work; /**< Work area for FFt's */ + dmat timedata_work; /**< Work area for timedata */ + vc j_vec_conj; /**< Work area for conjugate of j */ +} PowerSpectra; + +PowerSpectra* PowerSpectra_alloc(const us nfft, + const us nchannels, + const WindowType wt) { + + TRACE(15,"PowerSpectra_alloc"); + int rv; + + /* Check nfft */ + if(nfft % 2 != 0) { + WARN("nfft should be even"); + return NULL; + } + + /* ALlocate space */ + Fft* fft = Fft_alloc(nfft,nchannels); + 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); + ps->fft_work = cmat_alloc(nfft/2+1,nchannels); + ps->timedata_work= dmat_alloc(nfft,nchannels); + ps->j_vec_conj = vc_alloc(nfft/2+1); + + rv = window_create(wt,&ps->window,&ps->win_pow); + if(rv!=0) { + WARN("Error creating window function, continuing anyway"); + } + return ps; +} + +void PowerSpectra_free(PowerSpectra* ps) { + TRACE(15,"PowerSpectra_free"); + + Fft_free(ps->fft); + vd_free(&ps->window); + cmat_free(&ps->fft_work); + dmat_free(&ps->timedata_work); + vc_free(&ps->j_vec_conj); + a_free(ps); +} + + +int PowerSpectra_compute(const PowerSpectra* ps, + const dmat * timedata, + cmat * result) { + + TRACE(15,"PowerSpectra_compute"); + + const us nchannels = Fft_nchannels(ps->fft); + const us nfft = Fft_nfft(ps->fft); + uVARTRACE(15,nchannels); + const d win_pow = ps->win_pow; + dVARTRACE(15,win_pow); + + us i,j; + + /* Sanity checks for the matrices */ + dbgassert(timedata->n_cols == nchannels,"timedata n_cols " + "should be equal to nchannels"); + + 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 = ps->timedata_work; + + for(i=0;iwindow.data, + nfft); + + } + + /* print_dmat(&timedata_work); */ + + /* Compute fft of the time data */ + cmat fft_work = ps->fft_work; + Fft_fft(ps->fft, + &timedata_work, + &fft_work); + + + /* Scale fft such that power is easily comxputed */ + const c scale_fac = d_sqrt(2/win_pow)/nfft; + c_scale(fft_work.data,scale_fac,(nfft/2+1)*nchannels); + + for(i=0;i< nchannels;i++) { + /* Multiply DC term with 1/sqrt(2) */ + *getcmatval(&fft_work,0,i) *= 1/d_sqrt(2.)+0*I; + + /* Multiply Nyquist term with 1/sqrt(2) */ + *getcmatval(&fft_work,nfft/2,i) *= 1/d_sqrt(2.)+0*I; + } + + /* print_cmat(&fft_work); */ + + c* j_vec_conj = ps->j_vec_conj.data; + + /* Compute Cross-power spectra and store result */ + for(i =0; i nfft) { */ + +/* WARN("Overlap percentage results in offset of 0, or a too high number of */ +/* overlap" */ +/* " samples. Number of overlap samples should be < nfft"); */ + +/* WARN("Illegal overlap percentage. Should be 0 <= %% < 100"); */ +/* return NULL; */ +/* } */ + + +/* /\* ALlocate space *\/ */ +/* Fft fft; */ +/* rv = Fft_alloc(&fft,nfft,nchannels); */ +/* if(rv != SUCCESS) { */ +/* WARN("Fft allocation failed"); */ +/* return NULL; */ +/* } */ + +/* AvPowerSpectra* aps = a_malloc(sizeof(AvPowerSpectra)); */ +/* if(!aps) { */ +/* WARN("Allocation of AvPowerSpectra memory failed"); */ +/* return NULL; */ +/* } */ +/* ps->naverages = 0; */ +/* ps->overlap_offset = overlap_offset; */ + +/* /\* Allocate vectors and matrices *\/ */ +/* ps->prev_timedata = dmat_alloc(nfft,nchannels); */ + +/* return ps; */ +/* } */ +/* us AvPowerSpectra_getAverages(const AvPowerSpectra* ps) { */ +/* return ps->naverages; */ +/* } */ + +/* /\** */ +/* * Compute single power spectra by */ +/* * */ +/* * @param ps Initialized AvPowerSpectra structure */ +/* * @param timedata Timedata to compute for */ +/* * @param result Result */ +/* * */ +/* * @return */ +/* *\/ */ +/* static int AvPowerSpectra_computeSingle(const AvPowerSpectra* ps, */ +/* const dmat* timedata, */ +/* vc* results) { */ + +/* us nchannels = ps->fft.nchannels; */ +/* for(us channel=0;channelfft.nchannels; */ +/* const us nfft = ps->fft.nfft; */ + +/* dbgassert(timedata->n_cols == nchannels,"Invalid time data"); */ +/* dbgassert(timedata->n_rows == nfft,"Invalid time data"); */ + +/* dmat prevt = ps->prev_timedata; */ +/* us noverlap = ps->noverlap; */ + +/* if(ps->naverages != 0) { */ +/* /\* Copy the overlap buffer to the tbuf *\/ */ +/* copy_dmat_rows(&tbuf,&overlap,0,0,noverlap); */ + +/* /\* Copy the new timedata buffer to the tbuf *\/ */ +/* copy_dmat_rows(&tbuf,timedata,0,noverlap,(nfft-noverlap)); */ +/* } */ +/* else { */ +/* /\* Copy the overlap buffer to the tbuf *\/ */ +/* copy_dmat_rows(&tbuf,&overlap,0,0,noverlap); */ + + +/* } */ + +/* return SUCCESS; */ +/* } */ +/* void AvPowerSpectra_free(AvPowerSpectra* ps) { */ +/* TRACE(15,"AvPowerSpectra_free"); */ + +/* Fft_free(&ps->fft); */ +/* dmat_free(&ps->prev_timedata); */ +/* vd_free(&ps->window); */ + +/* a_free(ps); */ +/* } */ + +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/ps.h b/beamforming/c/ps.h new file mode 100644 index 0000000..9102dda --- /dev/null +++ b/beamforming/c/ps.h @@ -0,0 +1,67 @@ +// ps.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: Single sided power and cross-power spectra computation +// routines. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef PS_H +#define PS_H +#include "window.h" + +struct PowerSpectra_s; +struct AvPowerSpectra_s; + +typedef struct PowerSpectra_s PowerSpectra; +typedef struct AvPowerSpectra_s AvPowerSpectra; + + +PowerSpectra* PowerSpectra_alloc(const us nfft, + const us nchannels, + const WindowType wt); + +/** + * Compute power spectra, returns to a complex array + * + * @param[in] timedata Time data, should be of size nfft*nchannels + * @param [out] Cross-power spectra, array should be of size + * (nfft/2+1) x (nchannels*nchannels), such that the cross spectra of + * channel i with channel j at can be found as + * getvcval(0,i+j*nchannels). + * @return status code, SUCCESS on succes. + */ +int PowerSpectra_compute(const PowerSpectra*, + const dmat* timedata, + cmat *result); + +/** + * Free storage of the PowerSpectra + */ +void PowerSpectra_free(PowerSpectra*); + + +// AvPowerSpectra* AvPowerSpectra_alloc(const us nfft, +// const us nchannels, +// const d overlap_percentage, +// const WindowType wt); + + +/** + * Return the current number of averages taken to obtain the result. + * + * @return The number of averages taken. + */ +us AvPowerSpectra_getAverages(const AvPowerSpectra*); + + +int AvPowerSpectra_addTimeData(AvPowerSpectra*, + const dmat* timedata); + +/** + * Free storage of the AvPowerSpectra + */ +void AvPowerSpectra_free(AvPowerSpectra*); + +#endif // PS_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/signals.h b/beamforming/c/signals.h new file mode 100644 index 0000000..5ecf74a --- /dev/null +++ b/beamforming/c/signals.h @@ -0,0 +1,32 @@ +// signals.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// Several signal functions +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef SIGNALS_H +#define SIGNALS_H +#include "ascee_math.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;isize;i++) { + res+= d_pow(*getvdval(signal,i),2); + } + res /= signal->size; + return res; +} + + + +#endif // SIGNALS_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/tracer.c b/beamforming/c/tracer.c new file mode 100644 index 0000000..3d33850 --- /dev/null +++ b/beamforming/c/tracer.c @@ -0,0 +1,57 @@ +// si_tracer.c +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#if TRACER == 1 +#include +#include "tracer.h" +#include "types.h" + +#ifdef _REENTRANT +#include + +_Atomic(int) TRACERNAME = ATOMIC_VAR_INIT(DEFAULTTRACERLEVEL); + +void setTracerLevel(int level) { + atomic_store(&TRACERNAME,level); +} +int getTracerLevel() { + return atomic_load(&TRACERNAME); +} +#else + +int TRACERNAME; + +/* setTracerLevel and getTracerLevel are defined as macros in + * tracer.h */ +#endif + + + +void trace_impl(const char* file,int pos, const char * string){ + printf(annestr(TRACERNAME) ":%s:%i: %s\n",file,pos,string); +} +void fstrace_impl(const char* file,int pos,const char* fn){ + printf(annestr(TRACERNAME) ":%s:%i: start function: %s()\n",file,pos,fn); +} +void fetrace_impl(const char* file,int pos,const char* fn){ + printf(annestr(TRACERNAME) ":%s:%i: end function: %s()\n",file,pos,fn); +} +void ivartrace_impl(const char* pos,int line,const char* varname, int var){ + printf(annestr(TRACERNAME) ":%s:%i: %s = %i\n",pos,line,varname,var); +} +void uvartrace_impl(const char* pos,int line,const char* varname,size_t var){ + printf(annestr(TRACERNAME) ":%s:%i: %s = %zu\n",pos,line,varname,var); +} +void dvartrace_impl(const char* pos,int line,const char* varname, d var){ + printf(annestr(TRACERNAME) ":%s:%i: %s = %0.5e\n",pos,line,varname,var); +} +void cvartrace_impl(const char* pos,int line,const char* varname, c var){ + printf(annestr(TRACERNAME) ":%s:%i: %s = %0.5e+%0.5ei\n",pos,line,varname,creal(var),cimag(var)); +} +#endif + +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/tracer.h b/beamforming/c/tracer.h new file mode 100644 index 0000000..cd1fbc9 --- /dev/null +++ b/beamforming/c/tracer.h @@ -0,0 +1,219 @@ +// tracer.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// Basic tracing code for debugging. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef TRACER_H +#define TRACER_H + +// Some console colors +#define RESET "\033[0m" +#define BLACK "\033[30m" /* Black */ +#define RED "\033[31m" /* Red */ +#define GREEN "\033[32m" /* Green */ +#define YELLOW "\033[33m" /* Yellow */ +#define BLUE "\033[34m" /* Blue */ +#define MAGENTA "\033[35m" /* Magenta */ +#define CYAN "\033[36m" /* Cyan */ +#define WHITE "\033[37m" /* White */ +#define BOLDBLACK "\033[1m\033[30m" /* Bold Black */ +#define BOLDRED "\033[1m\033[31m" /* Bold Red */ +#define BOLDGREEN "\033[1m\033[32m" /* Bold Green */ +#define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */ +#define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */ +#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */ +#define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */ +#define BOLDWHITE "\033[1m\033[37m" /* Bold White */ + +#include +#include + +// Not so interesting part +#define rawstr(x) #x +#define namestr(x) rawstr(x) +#define annestr(x) namestr(x) +#define FILEWITHOUTPATH ( strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : strrchr(__FILE__, '\\') ? strrchr(__FILE__, '\\') + 1 : __FILE__ ) +// #define POS annestr(FILEWITHOUTPATH) ":" # __LINE__ << ": " +// End not so interesting part + +/** + * Produce a debug warning + */ +#define DBGWARN(a) \ + fprintf(stderr,RED);\ + fprintf(stderr,"%s(%d): ", \ + __FILE__, \ + __LINE__ \ + ); \ + fprintf(stderr,a); \ + fprintf(stderr,RESET "\n"); + +/** + * Produce a runtime warning + */ +#define WARN(a) \ + fprintf(stderr,RED);\ + fprintf(stderr,"WARNING: "); \ + fprintf(stderr,a); \ + fprintf(stderr,RESET "\n"); + +/** + * Fatal error, abort execution + */ +#define FATAL(a) \ + WARN(a); \ + abort(); + + + + +// **************************************** Tracer code +#ifndef TRACERPLUS +#define TRACERPLUS (0) +#endif + +// If PP variable TRACER is not defined, we automatically set it on. +#ifndef TRACER +#define TRACER 1 +#endif + +#if TRACER == 1 +#ifndef TRACERNAME + +#ifdef __GNUC__ +#warning TRACERNAME name not set, sol TRACERNAME set to 'defaulttracer' +#else +#pragma message("TRACERNAME name not set, sol TRACERNAME set to defaulttracer") +#endif + +#define TRACERNAME defaulttracer +#endif // ifndef TRACERNAME + +// Define this preprocessor definition to overwrite +// Use -O flag for compiler to remove the dead functions! +// In that case all cout's for TRACE() are removed from code +#ifndef DEFAULTTRACERLEVEL +#define DEFAULTTRACERLEVEL (15) +#endif + +#ifdef _REENTRANT +/** + * Set the tracer level at runtime + * + * @param level + */ +void setTracerLevel(int level); + +/** + * Obtain the tracer level + * + * @return level + */ +int getTracerLevel(); + +#else // Not reentrant +extern int TRACERNAME; +#define setTracerLevel(a) TRACERNAME = a; +static inline int getTracerLevel() { return TRACERNAME;} +#endif + +#include "types.h" + +// Use this preprocessor command to introduce one TRACERNAME integer per unit +/* Introduce one static logger */ +// We trust that the compiler will eliminate 'dead code', which means +// that if variable BUILDINTRACERLEVEL is set, the inner if statement +// will not be reached. +void trace_impl(const char* pos,int line,const char * string); +void fstrace_impl(const char* file,int pos,const char* fn); +void fetrace_impl(const char* file,int pos,const char* fn); +void dvartrace_impl(const char* pos,int line,const char* varname,d var); +void cvartrace_impl(const char* pos,int line,const char* varname,c var); +void ivartrace_impl(const char* pos,int line,const char* varname,int var); +void uvartrace_impl(const char* pos,int line,const char* varname,size_t var); + +/** + * Print a trace string + */ +#define TRACE(level,trace_string) \ + if (level+TRACERPLUS>=getTracerLevel()) \ + { \ + trace_impl(FILEWITHOUTPATH,__LINE__,trace_string ); \ + } + +/** + * Print start of function string + */ +#define fsTRACE(level) \ + if (level+TRACERPLUS>=getTracerLevel()) \ + { \ + fstrace_impl(FILEWITHOUTPATH,__LINE__, __FUNCTION__ ); \ + } + +/** + * Print end of function string + */ +#define feTRACE(level) \ + if (level+TRACERPLUS>=getTracerLevel()) \ + { \ + fetrace_impl(FILEWITHOUTPATH,__LINE__, __FUNCTION__ ); \ + } + + +/** + * Trace an int variable + */ +#define iVARTRACE(level,trace_var) \ + if (level+TRACERPLUS>=getTracerLevel()) \ + { \ + ivartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ + } + +/** + * Trace an unsigned int variable + */ +#define uVARTRACE(level,trace_var) \ + if (level+TRACERPLUS>=getTracerLevel()) \ + { \ + uvartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ + } +/** + * Trace a floating point value + */ +#define dVARTRACE(level,trace_var) \ + if (level+TRACERPLUS>=getTracerLevel()) \ + { \ + dvartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ + } +/** + * Trace a complex floating point value + */ +#define cVARTRACE(level,trace_var) \ + if (level+TRACERPLUS>=getTracerLevel()) \ + { \ + cvartrace_impl(FILEWITHOUTPATH,__LINE__,#trace_var,trace_var); \ + } + + +#else // TRACER !=1 +#define TRACE(l,a) +#define fsTRACE(l) +#define feTRACE(l) +#define setTracerLevel(a) +#define getTracerLevel() + +#define iVARTRACE(level,trace_var) +#define uVARTRACE(level,trace_var) +#define dVARTRACE(level,trace_var) +#define cVARTRACE(level,trace_var) + + +#endif // ######################################## TRACER ==1 + + + +#endif // TRACER_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/types.h b/beamforming/c/types.h new file mode 100644 index 0000000..042f8fb --- /dev/null +++ b/beamforming/c/types.h @@ -0,0 +1,43 @@ +// types.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// Typedefs and namespace pollution for stuff that is almost always +// needed. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef TYPES_H +#define TYPES_H + +/// We often use boolean values +#include // true, false +#include +#include +typedef size_t us; /* Size type I always use */ + +// To change the whole code to 32-bit floating points, change this to +// float. +#if ASCEE_FLOAT == 32 +typedef float d; /* Shortcut for double */ +typedef float complex c; +#elif ASCEE_FLOAT == 64 +typedef double d; /* Shortcut for double */ +typedef double complex c; +#else +#error ASCEE_FLOAT should be either 32 or 64 +#endif + + +/// I need these numbers so often, that they can be in the global +/// namespace. +#define SUCCESS 0 +#define INTERRUPTED (-3) +#define MALLOC_FAILED (-1) +#define FAILURE -2 + + + +#endif // ASCEE_TYPES_H +////////////////////////////////////////////////////////////////////// + diff --git a/beamforming/c/window.c b/beamforming/c/window.c new file mode 100644 index 0000000..3e2db16 --- /dev/null +++ b/beamforming/c/window.c @@ -0,0 +1,100 @@ +// window.c +// +// Author: J.A. de Jong -ASCEE +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "window.h" +#include "signals.h" + +/** + * Compute the Hann window + * + * @param i index + * @param N Number of indices + */ +static d hann(us i,us N) { + dbgassert(isize; + d (*win_fun)(us,us); + switch (wintype) { + case Hann: { + win_fun = hann; + break; + } + case Hamming: { + win_fun = hamming; + break; + } + case Blackman: { + win_fun = blackman; + break; + } + case Rectangular: { + win_fun = rectangle; + break; + } + default: + WARN("Unkown window function"); + return FAILURE; + break; + } + us index; + for(index=0;index +#include "ascee_assert.h" +#include "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[ASCEE_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 > ASCEE_MAX_NUM_THREADS) { + WARN("Number of workers too high in Workers_create"); + return NULL; + } + + dbgassert(init_fn,NULLPTRDEREF "init_fn"); + dbgassert(fn,NULLPTRDEREF "fn"); + dbgassert(free_fn,NULLPTRDEREF "free_fn"); + + Workers* w = a_malloc(sizeof(Workers)); + if(!w){ + WARN(ALLOCFAILED "Workers_create"); + return NULL; + } + + w->jq = jq; + w->w_alloc_fcn = init_fn; + w->fn = fn; + w->w_free_fcn = free_fn; + w->global_data = thread_global_data; + w->num_workers = num_workers; + + /* Initialize thread mutex */ + int rv = pthread_mutex_init(&w->global_data_mutex,NULL); + if(rv !=0) { + WARN("Mutex initialization failed"); + return NULL; + } + + + /* Create the threads */ + pthread_t* thread = w->worker_threads; + for(us i = 0; i < num_workers; i++) { + TRACE(15,"Creating thread"); + int rv = pthread_create(thread, + NULL, /* Thread attributes */ + threadfcn, /* Function */ + w); /* Data */ + if(rv!=0) { + WARN("Thread creation failed"); + return NULL; + } + thread++; + } + + return w; + +} + +void Workers_free(Workers* w) { + TRACE(15,"Workers_free"); + dbgassert(w,NULLPTRDEREF "w in Workers_free"); + dbgassert(w->jq,NULLPTRDEREF "w->jq in Workers_free"); + + for(us i=0;inum_workers;i++) { + /* Push the special NULL job. This will make the worker + * threads stop their execution. */ + JobQueue_push(w->jq,NULL); + } + + JobQueue_wait_alldone(w->jq); + +/* Join the threads */ + pthread_t* thread = w->worker_threads; + for(us i=0;inum_workers;i++) { + void* retval; + if(pthread_join(*thread,&retval)!=0) { + WARN("Error joining thread!"); + } + if((retval) != NULL) { + WARN("Thread returned with error status"); + } + thread++; + } + + /* Destroy the global data mutex */ + int rv = pthread_mutex_destroy(&w->global_data_mutex); + if(rv != 0){ + WARN("Mutex destroy failed. Do not know what to do."); + } + + /* All threads joined */ + a_free(w); + +} + +static void* threadfcn(void* thread_global_data) { + + TRACE(15,"Started worker thread function"); + dbgassert(thread_global_data,NULLPTRDEREF "thread_data in" + " threadfcn"); + Workers* w = (Workers*) thread_global_data; + + JobQueue* jq = w->jq; + worker_alloc_function walloc = w->w_alloc_fcn; + worker_free_function wfree = w->w_free_fcn; + worker_function worker_fn = w->fn; + void* global_data = w->global_data; + + dbgassert(jq,NULLPTRDEREF "jq in threadfcn"); + dbgassert(walloc,NULLPTRDEREF "walloc in threadfcn"); + dbgassert(wfree,NULLPTRDEREF "wfree in threadfcn"); + + int rv = pthread_mutex_lock(&w->global_data_mutex); + if(rv !=0) { + WARN("Global data mutex lock failed"); + pthread_exit((void*) 1); + } + + void* w_data = 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 */ + wfree(w_data); + TRACE(15,"Exiting thread. Goodbye"); + pthread_exit((void*) NULL); +} + + +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/c/worker.h b/beamforming/c/worker.h new file mode 100644 index 0000000..d16e215 --- /dev/null +++ b/beamforming/c/worker.h @@ -0,0 +1,61 @@ +// 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 WORKER_H +#define WORKER_H +#include "types.h" + +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 // WORKER_H +////////////////////////////////////////////////////////////////////// diff --git a/beamforming/cmake/FindCython.cmake b/beamforming/cmake/FindCython.cmake new file mode 100644 index 0000000..f1c4316 --- /dev/null +++ b/beamforming/cmake/FindCython.cmake @@ -0,0 +1,44 @@ +# Find the Cython compiler. +# +# This code sets the following variables: +# +# CYTHON_EXECUTABLE +# +# See also UseCython.cmake + +#============================================================================= +# Copyright 2011 Kitware, Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +#============================================================================= + +# Use the Cython executable that lives next to the Python executable +# if it is a local installation. +find_package( PythonInterp ) +if( PYTHONINTERP_FOUND ) + get_filename_component( _python_path ${PYTHON_EXECUTABLE} PATH ) + find_program( CYTHON_EXECUTABLE + NAMES cython cython.bat cython3 cython2 + HINTS ${_python_path} + ) +else() + find_program( CYTHON_EXECUTABLE + NAMES cython cython.bat cython3 + ) +endif() + + +include( FindPackageHandleStandardArgs ) +FIND_PACKAGE_HANDLE_STANDARD_ARGS( Cython REQUIRED_VARS CYTHON_EXECUTABLE ) + +mark_as_advanced( CYTHON_EXECUTABLE ) diff --git a/beamforming/cmake/ReplicatePythonSourceTree.cmake b/beamforming/cmake/ReplicatePythonSourceTree.cmake new file mode 100644 index 0000000..3d92d53 --- /dev/null +++ b/beamforming/cmake/ReplicatePythonSourceTree.cmake @@ -0,0 +1,6 @@ +# Note: when executed in the build dir, then CMAKE_CURRENT_SOURCE_DIR is the +# build dir. +file( COPY src/pyx/ DESTINATION "${CMAKE_ARGV3}" + FILES_MATCHING PATTERN "*.py" ) +file( COPY src/pyx/ DESTINATION "${CMAKE_ARGV3}" + FILES_MATCHING PATTERN "*.so" ) diff --git a/beamforming/cmake/UseCython.cmake b/beamforming/cmake/UseCython.cmake new file mode 100644 index 0000000..06fdf0c --- /dev/null +++ b/beamforming/cmake/UseCython.cmake @@ -0,0 +1,310 @@ +# Define a function to create Cython modules. +# +# For more information on the Cython project, see http://cython.org/. +# "Cython is a language that makes writing C extensions for the Python language +# as easy as Python itself." +# +# This file defines a CMake function to build a Cython Python module. +# To use it, first include this file. +# +# include( UseCython ) +# +# Then call cython_add_module to create a module. +# +# cython_add_module( ... ) +# +# To create a standalone executable, the function +# +# cython_add_standalone_executable( [MAIN_MODULE src1] ... ) +# +# To avoid dependence on Python, set the PYTHON_LIBRARY cache variable to point +# to a static library. If a MAIN_MODULE source is specified, +# the "if __name__ == '__main__':" from that module is used as the C main() method +# for the executable. If MAIN_MODULE, the source with the same basename as +# is assumed to be the MAIN_MODULE. +# +# Where is the name of the resulting Python module and +# ... are source files to be compiled into the module, e.g. *.pyx, +# *.py, *.c, *.cxx, etc. A CMake target is created with name . This can +# be used for target_link_libraries(), etc. +# +# The sample paths set with the CMake include_directories() command will be used +# for include directories to search for *.pxd when running the Cython complire. +# +# Cache variables that effect the behavior include: +# +# CYTHON_ANNOTATE +# CYTHON_NO_DOCSTRINGS +# CYTHON_FLAGS +# +# Source file properties that effect the build process are +# +# CYTHON_IS_CXX +# +# If this is set of a *.pyx file with CMake set_source_files_properties() +# command, the file will be compiled as a C++ file. +# +# See also FindCython.cmake + +#============================================================================= +# Copyright 2011 Kitware, Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +#============================================================================= + +# Configuration options. +set( CYTHON_ANNOTATE OFF + CACHE BOOL "Create an annotated .html file when compiling *.pyx." ) +set( CYTHON_NO_DOCSTRINGS OFF + CACHE BOOL "Strip docstrings from the compiled module." ) +set( CYTHON_FLAGS "" CACHE STRING + "Extra flags to the cython compiler." ) +mark_as_advanced( CYTHON_ANNOTATE CYTHON_NO_DOCSTRINGS CYTHON_FLAGS ) + +find_package( Cython REQUIRED ) + + +set( CYTHON_CXX_EXTENSION "cxx" ) +set( CYTHON_C_EXTENSION "c" ) + +# Create a *.c or *.cxx file from a *.pyx file. +# Input the generated file basename. The generate file will put into the variable +# placed in the "generated_file" argument. Finally all the *.py and *.pyx files. +function( compile_pyx _name generated_file ) + # Default to assuming all files are C. + set( cxx_arg "" ) + set( extension ${CYTHON_C_EXTENSION} ) + set( pyx_lang "C" ) + set( comment "Compiling Cython C source for ${_name}..." ) + + set( cython_include_directories "" ) + set( pxd_dependencies "" ) + set( pxi_dependencies "" ) + set( c_header_dependencies "" ) + set( pyx_locations "" ) + + foreach( pyx_file ${ARGN} ) + get_filename_component( pyx_file_basename "${pyx_file}" NAME_WE ) + + # Determine if it is a C or C++ file. + get_source_file_property( property_is_cxx ${pyx_file} CYTHON_IS_CXX ) + if( ${property_is_cxx} ) + set( cxx_arg "--cplus" ) + set( extension ${CYTHON_CXX_EXTENSION} ) + set( pyx_lang "CXX" ) + set( comment "Compiling Cython CXX source for ${_name}..." ) + endif() + + # Get the include directories. + get_source_file_property( pyx_location ${pyx_file} LOCATION ) + get_filename_component( pyx_path ${pyx_location} PATH ) + message(${pyx_path}) + get_directory_property( cmake_include_directories DIRECTORY ${pyx_path} INCLUDE_DIRECTORIES ) + list( APPEND cython_include_directories ${cmake_include_directories} ) + list( APPEND pyx_locations "${pyx_location}" ) + + # Determine dependencies. + # Add the pxd file will the same name as the given pyx file. + unset( corresponding_pxd_file CACHE ) + find_file( corresponding_pxd_file ${pyx_file_basename}.pxd + PATHS "${pyx_path}" ${cmake_include_directories} + NO_DEFAULT_PATH ) + if( corresponding_pxd_file ) + list( APPEND pxd_dependencies "${corresponding_pxd_file}" ) + endif() + + # pxd files to check for additional dependencies. + set( pxds_to_check "${pyx_file}" "${pxd_dependencies}" ) + set( pxds_checked "" ) + set( number_pxds_to_check 1 ) + while( ${number_pxds_to_check} GREATER 0 ) + foreach( pxd ${pxds_to_check} ) + list( APPEND pxds_checked "${pxd}" ) + list( REMOVE_ITEM pxds_to_check "${pxd}" ) + + # check for C header dependencies + file( STRINGS "${pxd}" extern_from_statements + REGEX "cdef[ ]+extern[ ]+from.*$" ) + foreach( statement ${extern_from_statements} ) + # Had trouble getting the quote in the regex + string( REGEX REPLACE "cdef[ ]+extern[ ]+from[ ]+[\"]([^\"]+)[\"].*" "\\1" header "${statement}" ) + unset( header_location CACHE ) + find_file( header_location ${header} PATHS ${cmake_include_directories} ) + if( header_location ) + list( FIND c_header_dependencies "${header_location}" header_idx ) + if( ${header_idx} LESS 0 ) + list( APPEND c_header_dependencies "${header_location}" ) + endif() + endif() + endforeach() + + # check for pxd dependencies + + # Look for cimport statements. + set( module_dependencies "" ) + file( STRINGS "${pxd}" cimport_statements REGEX cimport ) + foreach( statement ${cimport_statements} ) + if( ${statement} MATCHES from ) + string( REGEX REPLACE "from[ ]+([^ ]+).*" "\\1" module "${statement}" ) + else() + string( REGEX REPLACE "cimport[ ]+([^ ]+).*" "\\1" module "${statement}" ) + endif() + list( APPEND module_dependencies ${module} ) + endforeach() + list( REMOVE_DUPLICATES module_dependencies ) + # Add the module to the files to check, if appropriate. + foreach( module ${module_dependencies} ) + unset( pxd_location CACHE ) + find_file( pxd_location ${module}.pxd + PATHS "${pyx_path}" ${cmake_include_directories} NO_DEFAULT_PATH ) + if( pxd_location ) + list( FIND pxds_checked ${pxd_location} pxd_idx ) + if( ${pxd_idx} LESS 0 ) + list( FIND pxds_to_check ${pxd_location} pxd_idx ) + if( ${pxd_idx} LESS 0 ) + list( APPEND pxds_to_check ${pxd_location} ) + list( APPEND pxd_dependencies ${pxd_location} ) + endif() # if it is not already going to be checked + endif() # if it has not already been checked + endif() # if pxd file can be found + endforeach() # for each module dependency discovered + endforeach() # for each pxd file to check + list( LENGTH pxds_to_check number_pxds_to_check ) + endwhile() + + # Look for included pxi files + file(STRINGS "${pyx_file}" include_statements REGEX "include +['\"]([^'\"]+).*") + foreach(statement ${include_statements}) + string(REGEX REPLACE "include +['\"]([^'\"]+).*" "\\1" pxi_file "${statement}") + unset(pxi_location CACHE) + find_file(pxi_location ${pxi_file} + PATHS "${pyx_path}" ${cmake_include_directories} NO_DEFAULT_PATH) + if (pxi_location) + list(APPEND pxi_dependencies ${pxi_location}) + endif() + endforeach() # for each include statement found + + endforeach() # pyx_file + + # Set additional flags. + if( CYTHON_ANNOTATE ) + set( annotate_arg "--annotate" ) + endif() + + if( CYTHON_NO_DOCSTRINGS ) + set( no_docstrings_arg "--no-docstrings" ) + endif() + + if( "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "RelWithDebInfo" ) + set( cython_debug_arg "--gdb" ) + endif() + + if( "${PYTHONLIBS_VERSION_STRING}" MATCHES "^2." ) + set( version_arg "-2" ) + elseif( "${PYTHONLIBS_VERSION_STRING}" MATCHES "^3." ) + set( version_arg "-3" ) + else() + set( version_arg ) + endif() + + # Include directory arguments. + list( REMOVE_DUPLICATES cython_include_directories ) + set( include_directory_arg "" ) + foreach( _include_dir ${cython_include_directories} ) + set( include_directory_arg ${include_directory_arg} "-I" "${_include_dir}" ) + endforeach() + + # Determining generated file name. + set( _generated_file "${CMAKE_CURRENT_BINARY_DIR}/${_name}.${extension}" ) + set_source_files_properties( ${_generated_file} PROPERTIES GENERATED TRUE ) + set( ${generated_file} ${_generated_file} PARENT_SCOPE ) + + list( REMOVE_DUPLICATES pxd_dependencies ) + list( REMOVE_DUPLICATES c_header_dependencies ) + + # Add the command to run the compiler. + add_custom_command( OUTPUT ${_generated_file} + COMMAND ${CYTHON_EXECUTABLE} + ARGS ${cxx_arg} ${include_directory_arg} ${version_arg} + ${annotate_arg} ${no_docstrings_arg} ${cython_debug_arg} ${CYTHON_FLAGS} + --output-file ${_generated_file} ${pyx_locations} + DEPENDS ${pyx_locations} ${pxd_dependencies} ${pxi_dependencies} + IMPLICIT_DEPENDS ${pyx_lang} ${c_header_dependencies} + COMMENT ${comment} + ) + + # Remove their visibility to the user. + set( corresponding_pxd_file "" CACHE INTERNAL "" ) + set( header_location "" CACHE INTERNAL "" ) + set( pxd_location "" CACHE INTERNAL "" ) +endfunction() + +# cython_add_module( src1 src2 ... srcN ) +# Build the Cython Python module. +function( cython_add_module _name ) + set( pyx_module_sources "" ) + set( other_module_sources "" ) + foreach( _file ${ARGN} ) + if( ${_file} MATCHES ".*\\.py[x]?$" ) + list( APPEND pyx_module_sources ${_file} ) + else() + list( APPEND other_module_sources ${_file} ) + endif() + endforeach() + compile_pyx( ${_name} generated_file ${pyx_module_sources} ) + include_directories( ${PYTHON_INCLUDE_DIRS} ) + python_add_module( ${_name} ${generated_file} ${other_module_sources} ) + if( APPLE ) + set_target_properties( ${_name} PROPERTIES LINK_FLAGS "-undefined dynamic_lookup" ) + else() + target_link_libraries( ${_name} ${PYTHON_LIBRARIES} ) + endif() +endfunction() + +include( CMakeParseArguments ) +# cython_add_standalone_executable( _name [MAIN_MODULE src3.py] src1 src2 ... srcN ) +# Creates a standalone executable the given sources. +function( cython_add_standalone_executable _name ) + set( pyx_module_sources "" ) + set( other_module_sources "" ) + set( main_module "" ) + cmake_parse_arguments( cython_arguments "" "MAIN_MODULE" "" ${ARGN} ) + include_directories( ${PYTHON_INCLUDE_DIRS} ) + foreach( _file ${cython_arguments_UNPARSED_ARGUMENTS} ) + if( ${_file} MATCHES ".*\\.py[x]?$" ) + get_filename_component( _file_we ${_file} NAME_WE ) + if( "${_file_we}" STREQUAL "${_name}" ) + set( main_module "${_file}" ) + elseif( NOT "${_file}" STREQUAL "${cython_arguments_MAIN_MODULE}" ) + set( PYTHON_MODULE_${_file_we}_static_BUILD_SHARED OFF ) + compile_pyx( "${_file_we}_static" generated_file "${_file}" ) + list( APPEND pyx_module_sources "${generated_file}" ) + endif() + else() + list( APPEND other_module_sources ${_file} ) + endif() + endforeach() + + if( cython_arguments_MAIN_MODULE ) + set( main_module ${cython_arguments_MAIN_MODULE} ) + endif() + if( NOT main_module ) + message( FATAL_ERROR "main module not found." ) + endif() + get_filename_component( main_module_we "${main_module}" NAME_WE ) + set( CYTHON_FLAGS ${CYTHON_FLAGS} --embed ) + compile_pyx( "${main_module_we}_static" generated_file ${main_module} ) + add_executable( ${_name} ${generated_file} ${pyx_module_sources} ${other_module_sources} ) + target_link_libraries( ${_name} ${PYTHON_LIBRARIES} ${pyx_module_libs} ) +endfunction() diff --git a/beamforming/config.pxi.in b/beamforming/config.pxi.in new file mode 100644 index 0000000..4956f58 --- /dev/null +++ b/beamforming/config.pxi.in @@ -0,0 +1,42 @@ + +import numpy as np +cimport numpy as np +from libcpp cimport bool + +DEF ASCEE_FLOAT = "@ASCEE_FLOAT@" + +IF ASCEE_FLOAT == "double": + ctypedef double d + ctypedef double complex c + NUMPY_FLOAT_TYPE = np.float64 + NUMPY_COMPLEX_TYPE = np.complex128 +ELSE: + ctypedef float d + ctypedef float complex c + NUMPY_FLOAT_TYPE = np.float32 + NUMPY_COMPLEX_TYPE = np.complex128 + +ctypedef size_t us + +cdef extern from "math.h": + ctypedef struct dmat: + d* data + us n_cols,n_rows + ctypedef struct cmat: + c* data + us n_cols,n_rows + +# cdef np.ndarray empty(us nrows,us ncols,dtype): +cdef dmat dmat_from_array(d[::1,:] arr): + cdef dmat result + result.data = &arr[0,0] + result.n_rows = arr.shape[0] + result.n_cols = arr.shape[1] + return result + +cdef cmat cmat_from_array(c[::1,:] arr): + cdef cmat result + result.data = &arr[0,0] + result.n_rows = arr.shape[0] + result.n_cols = arr.shape[1] + return result diff --git a/beamforming/fft.pyx b/beamforming/fft.pyx new file mode 100644 index 0000000..f4ef18d --- /dev/null +++ b/beamforming/fft.pyx @@ -0,0 +1,95 @@ +include "config.pxi" + +cdef extern from "fft.h": + ctypedef struct c_Fft "Fft" + c_Fft* Fft_alloc(us nfft,us nchannels) + void Fft_free(c_Fft*) + void Fft_fft(c_Fft*,dmat * timedate,cmat * res) nogil + us Fft_nchannels(c_Fft*) + us Fft_nfft(c_Fft*) + +cdef class Fft: + cdef: + c_Fft* _fft + + def __cinit__(self, us nfft,us nchannels): + self._fft = Fft_alloc(nfft,nchannels) + if self._fft == NULL: + raise RuntimeError('Fft allocation failed') + + def __dealloc__(self): + if self._fft!=NULL: + Fft_free(self._fft) + + def fft(self,d[::1,:] timedata): + + cdef us nfft = Fft_nfft(self._fft) + cdef us nchannels = Fft_nchannels(self._fft) + assert timedata.shape[0] ==nfft + assert timedata.shape[1] == nchannels + + result = np.empty((nfft//2+1, + nchannels), + dtype=NUMPY_COMPLEX_TYPE,order='F') + # result[:,:] = np.nan+1j*np.nan + + cdef cmat r = cmat_from_array(result) + cdef dmat t = dmat_from_array(timedata) + Fft_fft(self._fft,&t,&r) + + return result + + +cdef extern from "window.h": + ctypedef enum WindowType: + Hann = 0 + Hamming = 1 + Blackman = 2 + Rectangular = 3 + +cdef extern from "ps.h": + ctypedef struct c_PowerSpectra "PowerSpectra" + c_PowerSpectra* PowerSpectra_alloc(const us nfft, + const us nchannels, + const WindowType wt) + int PowerSpectra_compute(const c_PowerSpectra* ps, + const dmat * timedata, + cmat * result) + + + void PowerSpectra_free(c_PowerSpectra*) + +# cdef class PowerSpectra: +# cdef: +# c_PowerSpectra* _ps + +# def __cinit__(self, us nfft,us nchannels): +# self._ps = PowerSpectra_alloc(nfft,nchannels,Rectangular) +# if self._ps == NULL: +# raise RuntimeError('PowerSpectra allocation failed') + +# def compute(self,d[::1,:] timedata): +# cdef: +# us nchannels = timedata.shape[1] +# us nfft = timedata.shape[0] +# int rv +# dmat td +# cmat result_mat +# td.data = &timedata[0,0] +# td.n_rows = nfft +# td.n_cols = nchannels +# result = np.empty((nfft//2+1,nchannels*nchannels), +# dtype = NUMPY_COMPLEX_TYPE, +# order='F') +# result_mat = cmat_from_array(result) + +# rv = PowerSpectra_compute(self._ps,&td,&result_mat) +# if rv !=0: +# raise RuntimeError("Error computing power spectra") + +# return result + + +# def __dealloc__(self): +# if self._ps != NULL: +# PowerSpectra_free(self._ps) diff --git a/fftpack/CMakeLists.txt b/fftpack/CMakeLists.txt new file mode 100644 index 0000000..5c4f22b --- /dev/null +++ b/fftpack/CMakeLists.txt @@ -0,0 +1,8 @@ +set(fftpack_src numpy/fftpack.c) + +add_library(fftpack + fftpack.c + ) + + +target_link_libraries(fftpack m) diff --git a/fftpack/fftpack.c b/fftpack/fftpack.c new file mode 100644 index 0000000..eaf172c --- /dev/null +++ b/fftpack/fftpack.c @@ -0,0 +1,1501 @@ +/* + * fftpack.c : A set of FFT routines in C. + * Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985). +*/ +#define NPY_VISIBILITY_HIDDEN +/* #define NPY_NO_DEPRECATED_API NPY_API_VERSION */ + +#include +#include + + +#define DOUBLE +#ifdef DOUBLE +#define Treal double +#else +#define Treal float +#endif + + +#define ref(u,a) u[a] + +#define MAXFAC 13 /* maximum number of factors in factorization of n */ +#define NSPECIAL 4 /* number of factors for which we have special-case routines */ + +#ifdef __cplusplus +extern "C" { +#endif + + +/* ---------------------------------------------------------------------- + passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd. +----------------------------------------------------------------------- */ + +static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign) + /* isign==+1 for backward transform */ + { + int i, k, ah, ac; + Treal ti2, tr2; + if (ido <= 2) { + for (k=0; k= l1) { + for (j=1; j idp) idlj -= idp; + war = wa[idlj - 2]; + wai = wa[idlj-1]; + for (ik=0; ik= l1) { + for (j=1; j= l1) { + for (k=0; k= l1) { + for (j=1; j= l1) { + for (k=0; k= l1) { + for (j=1; j= l1) { + for (j=1; j 5) { + wa[i1-1] = wa[i-1]; + wa[i1] = wa[i]; + } + } + l1 = l2; + } + } /* cffti1 */ + + +NPY_VISIBILITY_HIDDEN void npy_cffti(int n, Treal wsave[]) + { + int iw1, iw2; + if (n == 1) return; + iw1 = 2*n; + iw2 = iw1 + 2*n; + cffti1(n, wsave+iw1, (int*)(wsave+iw2)); + } /* npy_cffti */ + + /* ------------------------------------------------------------------- +rfftf1, rfftb1, npy_rfftf, npy_rfftb, rffti1, npy_rffti. Treal FFTs. +---------------------------------------------------------------------- */ + +static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2]) + { + int i; + int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1; + Treal *cinput, *coutput; + nf = ifac[1]; + na = 1; + l2 = n; + iw = n-1; + for (k1 = 1; k1 <= nf; ++k1) { + kh = nf - k1; + ip = ifac[kh + 2]; + l1 = l2 / ip; + ido = n / l2; + idl1 = ido*l1; + iw -= (ip - 1)*ido; + na = !na; + if (na) { + cinput = ch; + coutput = c; + } else { + cinput = c; + coutput = ch; + } + switch (ip) { + case 4: + ix2 = iw + ido; + ix3 = ix2 + ido; + radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]); + break; + case 2: + radf2(ido, l1, cinput, coutput, &wa[iw]); + break; + case 3: + ix2 = iw + ido; + radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]); + break; + case 5: + ix2 = iw + ido; + ix3 = ix2 + ido; + ix4 = ix3 + ido; + radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); + break; + default: + if (ido == 1) + na = !na; + if (na == 0) { + radfg(ido, ip, l1, idl1, c, ch, &wa[iw]); + na = 1; + } else { + radfg(ido, ip, l1, idl1, ch, c, &wa[iw]); + na = 0; + } + } + l2 = l1; + } + if (na == 1) return; + for (i = 0; i < n; i++) c[i] = ch[i]; + } /* rfftf1 */ + + +static void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2]) + { + int i; + int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; + Treal *cinput, *coutput; + nf = ifac[1]; + na = 0; + l1 = 1; + iw = 0; + for (k1=1; k1<=nf; k1++) { + ip = ifac[k1 + 1]; + l2 = ip*l1; + ido = n / l2; + idl1 = ido*l1; + if (na) { + cinput = ch; + coutput = c; + } else { + cinput = c; + coutput = ch; + } + switch (ip) { + case 4: + ix2 = iw + ido; + ix3 = ix2 + ido; + radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]); + na = !na; + break; + case 2: + radb2(ido, l1, cinput, coutput, &wa[iw]); + na = !na; + break; + case 3: + ix2 = iw + ido; + radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]); + na = !na; + break; + case 5: + ix2 = iw + ido; + ix3 = ix2 + ido; + ix4 = ix3 + ido; + radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); + na = !na; + break; + default: + radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]); + if (ido == 1) na = !na; + } + l1 = l2; + iw += (ip - 1)*ido; + } + if (na == 0) return; + for (i=0; i0,"Invalid nfft") + dbgassert(fftr,ALLOCFAILED "Fftr_alloc"); + fftr->work_area = a_malloc(sizeof(d)*(3*nfft+15)); + fftr->wsave_ptr = &fftr->work_area[nfft]; + fftr->nfft = nfft; + #if ASCEE_DOUBLE_PRECISION + // dffti_(&nfft,fftr->wsave); + npy_rffti(nfft,fftr->wsave_ptr); + #else + // rffti_(&nfft,fftr->wsave); + #endif + + feTRACE(15); + return fftr; +} +void Fftr_free(Fftr* fftr) { + dbgassert(fftr,NULLPTRDEREF "Fftr_free"); + a_free(fftr->work_area); + a_free(fftr); +} + + +/** + * + * + * @param nfft + * @param wsave + * @param input + * @param work + * @param result + */ +static inline void Fftr_fftr(Fftr* fftr,d* input,c* result) { + + + // Copy contents of input to the work area + d_copy(fftr->work_area,input,fftr->nfft); + + int nfft = fftr->nfft; + + #if ASCEE_DOUBLE_PRECISION + // dfftf_(&nfft,fftr->work,fftr->wsave); + npy_rfftf(nfft,fftr->work_area,fftr->wsave_ptr); + #else + NOT TESTED + rfftf_(&nfft,fftr->work_area,fftr->wsave_ptr); + #endif + + result[0] = fftr->work_area[0]; + d_copy((d*) (&result[1]),&fftr->work_area[1],nfft); + // Not portable way of setting imaginary part to zero. Works with + // gcc, though. + __imag__ result[nfft/2] = 0; + +} + + + +#endif // FFTPACK_H +////////////////////////////////////////////////////////////////////// diff --git a/fftpack/npy_fftpack.h b/fftpack/npy_fftpack.h new file mode 100644 index 0000000..aa706ea --- /dev/null +++ b/fftpack/npy_fftpack.h @@ -0,0 +1,29 @@ +/* + * This file is part of tela the Tensor Language. + * Copyright (c) 1994-1995 Pekka Janhunen + */ + +#ifdef __cplusplus +extern "C" { +#endif +#define NPY_VISIBILITY_HIDDEN + +#define DOUBLE + +#ifdef DOUBLE +#define Treal double +#else +#define Treal float +#endif + +extern NPY_VISIBILITY_HIDDEN void npy_cfftf(int N, Treal data[], const Treal wrk[]); +extern NPY_VISIBILITY_HIDDEN void npy_cfftb(int N, Treal data[], const Treal wrk[]); +extern NPY_VISIBILITY_HIDDEN void npy_cffti(int N, Treal wrk[]); + +extern NPY_VISIBILITY_HIDDEN void npy_rfftf(int N, Treal data[], const Treal wrk[]); +extern NPY_VISIBILITY_HIDDEN void npy_rfftb(int N, Treal data[], const Treal wrk[]); +extern NPY_VISIBILITY_HIDDEN void npy_rffti(int N, Treal wrk[]); + +#ifdef __cplusplus +} +#endif diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..0b4c640 --- /dev/null +++ b/setup.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +@author: J.A. de Jong - ASCEE +""" +from setuptools import setup + +descr = "Python wrappers around several C++ optimized beamforming codes" + + +setup( + name="Beamforming", + version="1.0", + packages=['beamforming'], + author='J.A. de Jong - ASCEE', + author_email="j.a.dejong@ascee.nl", + # Project uses reStructuredText, so ensure that the docutils get + # installed or upgraded on the target machine + install_requires=['matplotlib>=1.0', 'scipy>=0.17', 'numpy>=1.0'], + license='MIT', + + description=descr, + keywords="Beamforming", + url="http://www.ascee.nl/beamforming/", # project home page, if any + +) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt new file mode 100644 index 0000000..79afbda --- /dev/null +++ b/test/CMakeLists.txt @@ -0,0 +1,7 @@ +include_directories(${CMAKE_SOURCE_DIR}/beamforming/c) +add_executable(test_bf test_bf.c) +add_executable(test_workers test_workers.c) +add_executable(test_fft test_fft.c) +target_link_libraries(test_bf beamforming_lib) +target_link_libraries(test_fft beamforming_lib pthread) +target_link_libraries(test_workers beamforming_lib pthread) diff --git a/test/benchmark.py b/test/benchmark.py new file mode 100644 index 0000000..1039303 --- /dev/null +++ b/test/benchmark.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Jan 15 19:45:33 2018 + +@author: anne +""" +import numpy as np +from beamforming.fft import Fft + +nfft=2**17 + +print('nfft:',nfft) +nchannels = 50 +number_run = 10 + +t = np.linspace(0,1,nfft+1) + +# Using transpose to get the strides right +x = np.random.randn(nchannels,nfft).T + + +import time +start = time.time() +for i in range(number_run): + X = np.fft.rfft(x,axis=0) +end = time.time() +print("Time numpy fft:",end-start) +# X =np.fft.fft(x) +#X =np.fft.rfft(x) + +fft = Fft(nfft,nchannels) + +start = time.time() +for i in range(number_run): + # print('--- run %i' %i) + fft.fft(x) +end = time.time() + +print("Time ASCEE fft:",end-start) + + diff --git a/test/fft_test.py b/test/fft_test.py new file mode 100644 index 0000000..e2a0012 --- /dev/null +++ b/test/fft_test.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Jan 15 19:45:33 2018 + +@author: anne +""" +import numpy as np +from beamforming.fft import Fft + +nfft=6 +print('nfft:',nfft) +print(nfft) +nchannels = 1 + +t = np.linspace(0,1,nfft+1)[:-1] +# print(t) +x1 = 1+np.sin(2*np.pi*t)+3.2*np.cos(2*np.pi*t)+np.sin(7*np.pi*t) +x = np.vstack([x1.T]*nchannels).T +# Using transpose to get the strides right +x = np.random.randn(nchannels,nfft).T +# x.strides = (8,nfft*8)x +# print("signal:",x) + +X = np.fft.rfft(x,axis=0) +print('Numpy fft') +print(X) + +fft = Fft(nfft,nchannels) +Y = fft.fft(x) +print('Fftpack fft') +print(Y) +print('end python script') diff --git a/test/ps_test.py b/test/ps_test.py new file mode 100644 index 0000000..8acb94e --- /dev/null +++ b/test/ps_test.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Jan 15 19:45:33 2018 + +@author: anne +""" +import numpy as np +from beamforming.fft import Fft, PowerSpectra + +nfft=2**10 +print('nfft:',nfft) +print(nfft) +nchannels = 8 + +t = np.linspace(0,1,nfft+1) +# print(t) +x1 = (np.cos(4*np.pi*t[:-1])+3.2*np.sin(6*np.pi*t[:-1]))[:,np.newaxis]+10 +x = np.vstack([x1.T]*nchannels).T +# Using transpose to get the strides right +x = np.random.randn(nchannels,nfft).T +print("strides: ",x.strides) +# x.strides = (8,nfft*8)x +# print("signal:",x) + +xms = np.sum(x**2,axis=0)/nfft +print('Total signal power time domain: ', xms) + +X = np.fft.rfft(x,axis=0) +# X =np.fft.fft(x) +#X =np.fft.rfft(x) + +# print(X) +Xs = 2*X/nfft +Xs[np.where(np.abs(Xs) < 1e-10)] = 0 +Xs[0] /= np.sqrt(2) +Xs[-1] /= np.sqrt(2) +# print('single sided amplitude spectrum:\n',Xs) + + +power = Xs*np.conj(Xs)/2 + +# print('Frequency domain signal power\n', power) +print('Total signal power', np.sum(power,axis=0).real) + +pstest = PowerSpectra(nfft,nchannels) +ps = pstest.compute(x) + +fft = Fft(nfft,nchannels) +fft.fft(x) + +ps[np.where(np.abs(ps) < 1e-10)] = 0+0j + +print('our ps: \n' , ps) + +print('Our total signal power: ',np.sum(ps,axis=0).real) diff --git a/test/test_bf.c b/test/test_bf.c new file mode 100644 index 0000000..98779d5 --- /dev/null +++ b/test/test_bf.c @@ -0,0 +1,46 @@ +// test_bf.c +// +// Author: J.A. de Jong -ASCEE +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#include "ascee_math.h" + +int main() { + + iVARTRACE(15,getTracerLevel()); + /* vd vec1 = vd_alloc(3); */ + /* vd_set(&vec1,2); */ + + /* vd vec2 = vd_alloc(3); */ + /* vd_set(&vec2,3); */ + + /* print_vd(&vec1); */ + + /* vd res = vd_alloc(3); */ + /* d_elem_prod_d(res.data,vec1.data,vec2.data,3); */ + + /* print_vd(&res); */ + + + vc vc1 = vc_alloc(3); + vc_set(&vc1,2+2I); + print_vc(&vc1); + + vc vc2 = vc_alloc(3); + vc_set(&vc2,2-2I); + setvecval(&vc2,0,10); + print_vc(&vc2); + + + vc res2 = vc_alloc(3); + c_elem_prod_c(res2.data,vc1.data,vc2.data,3); + + print_vc(&res2); + + + +} + +////////////////////////////////////////////////////////////////////// diff --git a/test/test_fft.c b/test/test_fft.c new file mode 100644 index 0000000..5df900d --- /dev/null +++ b/test/test_fft.c @@ -0,0 +1,31 @@ +// test_bf.c +// +// Author: J.A. de Jong -ASCEE +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "fft.h" +#include "tracer.h" + + +int main() { + + iVARTRACE(15,getTracerLevel()); + + Fft* fft = Fft_alloc(100000,8); + + Fft_fft(fft,NULL,NULL); + + Fft_free(fft); + + return 0; +} + + + + +////////////////////////////////////////////////////////////////////// + + diff --git a/test/test_workers.c b/test/test_workers.c new file mode 100644 index 0000000..321412d --- /dev/null +++ b/test/test_workers.c @@ -0,0 +1,70 @@ +// test_bf.c +// +// Author: J.A. de Jong -ASCEE +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#include "worker.h" +#include "mq.h" +#include "tracer.h" + +static void* walloc(void*); +static int worker(void*,void*); +static void wfree(void*); + + +int main() { + + fsTRACE(15); + + iVARTRACE(15,getTracerLevel()); + + int njobs = 4; + JobQueue* jq = JobQueue_alloc(njobs); + assert(jq); + + Workers* w = Workers_create(njobs, + jq, + walloc, + worker, + wfree, + (void*) 101); + assert(w); + + for(us i=0; i< njobs; i++) { + iVARTRACE(15,i); + JobQueue_push(jq,(void*) i+1); + } + + JobQueue_wait_alldone(jq); + Workers_free(w); + JobQueue_free(jq); + + return 0; +} + +static void* walloc(void* data) { + TRACE(15,"WALLOC"); + uVARTRACE(15,(int) data); + return (void*) getpid(); +} + +static int worker(void* w_data,void* tj) { + + TRACE(15,"worker"); + + sleep(4); + + return 0; + +} +static void wfree(void* w_data) { + TRACE(15,"wfree"); +} + + + +////////////////////////////////////////////////////////////////////// + +