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

Renamed lasp_math.h to lasp_mat.h as it defines matrices, not math. Made a typedef from vd to dmat such that the types of these are the same. Same goes for vc and cmat

Browse Source
This commit is contained in:
Anne de Jong 2018-04-01 14:01:07 +02:00 committed by Anne de Jong
parent c8332f7f74
commit ce7f31ca86
23 changed files with 268 additions and 313 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
lasp/c/CMakeLists.txt
View File

@ -5,7 +5,7 @@ endif(!LASP_DEBUG)
add_library(lasp_lib
lasp_fft.c
lasp_math.c
lasp_mat.c
lasp_math_raw.c
lasp_alg.c
lasp_assert.c

8
lasp/c/lasp_alg.c
View File

@ -8,10 +8,12 @@
#include "lasp_alg.h"
void cmv_dot(const cmat* A,const vc* restrict x,vc* restrict b){
dbgassert(A && x && b,NULLPTRDEREF);
assert_vx(b);
assert_vx(x);
dbgassert(A->n_cols == x->n_rows,SIZEINEQUAL);
dbgassert(A->n_rows == b->n_rows,SIZEINEQUAL);
dbgassert(A->n_rows == b->size,SIZEINEQUAL);
dbgassert(A->n_cols == x->size,SIZEINEQUAL);
#if LASP_USE_BLAS == 1
dbgassert(false,"Untested function. Is not functional for strides");
/* typedef enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER; */

20
lasp/c/lasp_alg.h
View File

@ -8,7 +8,7 @@
#pragma once
#ifndef LASP_ALG_H
#define LASP_ALG_H
#include "lasp_math.h"
#include "lasp_mat.h"
/**
* Compute the dot product of two vectors of floats
@ -19,8 +19,10 @@
*/
static inline d vd_dot(const vd * a,const vd* b) {
dbgassert(a && b,NULLPTRDEREF);
dbgassert(a->size == b->size,SIZEINEQUAL);
return d_dot(getvdval(a,0),getvdval(b,0),a->size);
assert_vx(a);
assert_vx(b);
dbgassert(a->n_rows == b->n_rows,SIZEINEQUAL);
return d_dot(getvdval(a,0),getvdval(b,0),a->n_rows);
}
/**
@ -93,11 +95,10 @@ static inline void cmat_add_cmat(cmat* x,cmat* y,c fac) {
*/
static inline void vd_elem_prod(vd* result,const vd* a,const vd* b) {
dbgassert(result && a && b,NULLPTRDEREF);
dbgassert(result->size==a->size,SIZEINEQUAL);
dbgassert(b->size==a->size,SIZEINEQUAL);
assert_equalsize(a,b);
d_elem_prod_d(getvdval(result,0),
getvdval(a,0),
getvdval(b,0),a->size);
getvdval(b,0),a->n_rows);
}
/**
* Compute the element-wise (Hadamard) product of a and b, and store
@ -110,12 +111,13 @@ static inline void vd_elem_prod(vd* result,const vd* a,const vd* b) {
static inline void vc_hadamard(vc* result,const vc* a,const vc* b) {
fsTRACE(15);
dbgassert(result && a && b,NULLPTRDEREF);
dbgassert(result->size==a->size,SIZEINEQUAL);
dbgassert(b->size==a->size,SIZEINEQUAL);
assert_equalsize(a,b);
assert_equalsize(a,result);
c_hadamard(getvcval(result,0),
getvcval(a,0),
getvcval(b,0),
a->size);
a->n_rows);
check_overflow_vx(*result);
check_overflow_vx(*a);

10
lasp/c/lasp_aps.c
View File

@ -64,7 +64,7 @@ void AvPowerSpectra_free(AvPowerSpectra* aps) {
dFifo_free(aps->fifo);
dmat_free(&aps->buffer);
us nweight = aps->weighting.size;
us nweight = aps->weighting.n_rows;
if(nweight > 0) {
for(us blockno = 0; blockno < nweight; blockno++) {
cmat_free(&aps->ps_storage[blockno]);
@ -127,7 +127,7 @@ AvPowerSpectra* AvPowerSpectra_alloc(const us nfft,
aps->oldest_block = 0;
if(weighting) {
us nweight = weighting->size;
us nweight = weighting->n_rows;
iVARTRACE(15,nweight);
/* Allocate vectors and matrices */
aps->ps_storage = a_malloc(nweight*sizeof(cmat));
@ -137,12 +137,12 @@ AvPowerSpectra* AvPowerSpectra_alloc(const us nfft,
}
/* Allocate space and copy weighting coefficients */
aps->weighting = vd_alloc(weighting->size);
aps->weighting = vd_alloc(weighting->n_rows);
vd_copy(&aps->weighting,weighting);
}
else {
TRACE(15,"no weighting");
aps->weighting.size = 0;
aps->weighting.n_rows = 0;
}
aps->ps_result = cmat_alloc(nfft/2+1,nchannels*nchannels);
@ -186,7 +186,7 @@ static void AvPowerSpectra_addBlock(AvPowerSpectra* aps,
ps_single);
vd weighting = aps->weighting;
us nweight = weighting.size;
us nweight = weighting.n_rows;
if(nweight == 0) {

2
lasp/c/lasp_aps.h
View File

@ -9,7 +9,7 @@
#ifndef LASP_APS_H
#define LASP_APS_H
#include "lasp_types.h"
#include "lasp_math.h"
#include "lasp_mat.h"
#include "lasp_window.h"
typedef enum {

2
lasp/c/lasp_decimation.c
View File

@ -62,7 +62,7 @@ Decimator* Decimator_create(us nchannels,DEC_FAC df) {
dec->nchannels = nchannels;
dec->dec_fac = filter->dec_fac;
dmat h = dmat_foreign(filter->ntaps,1,filter->h);
dmat h = dmat_foreign_data(filter->ntaps,1,filter->h,false);
for(us channelno=0;channelno<nchannels;channelno++) {
dec->fbs[channelno] = FilterBank_create(&h,DEC_FFT_LEN);

2
lasp/c/lasp_decimation.h
View File

@ -10,7 +10,7 @@
#ifndef LASP_DECIMATION_H
#define LASP_DECIMATION_H
#include "lasp_types.h"
#include "lasp_math.h"
#include "lasp_mat.h"
typedef struct Decimator_s Decimator;

2
lasp/c/lasp_dfifo.h
View File

@ -9,7 +9,7 @@
#ifndef LASP_DFIFO_H
#define LASPDFIFO_H
#include "lasp_types.h"
#include "lasp_math.h"
#include "lasp_mat.h"
typedef struct dFifo_s dFifo;

10
lasp/c/lasp_fft.c
View File

@ -54,11 +54,11 @@ void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) {
fsTRACE(15);
dbgassert(fft && freqdata && result,NULLPTRDEREF);
const us nfft = fft->nfft;
dbgassert(result->size == nfft,
dbgassert(result->n_rows == nfft,
"Invalid size for time data rows."
" Should be equal to nfft");
dbgassert(freqdata->size == (nfft/2+1),"Invalid number of rows in"
dbgassert(freqdata->n_rows == (nfft/2+1),"Invalid number of rows in"
" result array");
@ -113,11 +113,13 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
dbgassert(fft && timedata && result,NULLPTRDEREF);
const us nfft = fft->nfft;
dbgassert(timedata->size == nfft,
assert_vx(timedata);
assert_vx(result);
dbgassert(timedata->n_rows == nfft,
"Invalid size for time data rows."
" Should be equal to nfft");
dbgassert(result->size == (nfft/2+1),"Invalid number of rows in"
dbgassert(result->n_rows == (nfft/2+1),"Invalid number of rows in"
" result array");
d* result_ptr = (d*) getvcval(result,0);

2
lasp/c/lasp_fft.h
View File

@ -9,7 +9,7 @@
#ifndef LASP_FFT_H
#define LASP_FFT_H
#include "lasp_types.h"
#include "lasp_math.h"
#include "lasp_mat.h"
/**
* Perform forward FFT's on real time data.

4
lasp/c/lasp_filterbank.c
View File

@ -100,9 +100,7 @@ dmat FilterBank_filter(FilterBank* fb,
const us nfilters = fb->filters.n_cols;
/* Push samples to the input fifo */
dmat samples = dmat_foreign_vd((vd*) x);
dFifo_push(fb->input_fifo,&samples);
dmat_free(&samples);
dFifo_push(fb->input_fifo,x);
dmat input_block = dmat_alloc(nfft,1);

2
lasp/c/lasp_filterbank.h
View File

@ -13,7 +13,7 @@
#ifndef LASP_FILTERBANK_H
#define LASP_FILTERBANK_H
#include "lasp_types.h"
#include "lasp_math.h"
#include "lasp_mat.h"
typedef struct FilterBank_s FilterBank;
/**

35
lasp/c/lasp_math.c → lasp/c/lasp_mat.c
View File

@ -1,4 +1,4 @@
// lasp_math.c
// lasp_mat.c
//
// Author: J.A. de Jong -ASCEE
//
@ -6,10 +6,8 @@
//
//////////////////////////////////////////////////////////////////////
#define TRACERPLUS (-10)
#include "lasp_math.h"
#include "lasp_mat.h"
#include "lasp_assert.h"
#include "lasp_math.h"
#include "lasp_tracer.h"
#include <math.h>
@ -50,35 +48,6 @@ void print_cmat(const cmat* m) {
}
feTRACE(50);
}
void print_vc(const vc* m) {
fsTRACE(50);
size_t row;
for(row=0;row<m->size;row++){
c val = *getvcval(m,row);
d rval = creal(val);
d ival = cimag(val);
indent_trace();
printf("%c%2.2e%c%2.2ei ",rval< 0 ?'-': ' ', d_abs(rval),ival<0 ? '-' : '+',d_abs(ival) ) ;
printf("\n");
}
feTRACE(50);
}
void print_vd(const vd* m) {
fsTRACE(50);
size_t row;
iVARTRACE(20,m->size);
for(row=0;row<m->size;row++){
d rval = *getvdval(m,row);
indent_trace();
printf("%c%2.2e ",rval< 0 ? '\r': ' ',rval);
printf("\n");
}
feTRACE(50);
}
#endif

358
lasp/c/lasp_math.h → lasp/c/lasp_mat.h
View File

@ -1,4 +1,4 @@
// lasp_math.h
// lasp_mat.h
//
// Author: J.A. de Jong - ASCEE
//
@ -6,29 +6,13 @@
// copying of matrices and vectors.
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_MATH_H
#define LASP_MATH_H
#ifndef LASP_MAT_H
#define LASP_MAT_H
#include "lasp_math_raw.h"
#include "lasp_alloc.h"
#include "lasp_tracer.h"
#include "lasp_assert.h"
/// Vector of floating point numbers
typedef struct {
us size;
bool _foreign_data;
d* _data; /**< Pointer set if data storage is
intern. If this is set to zero, the
vector is a sub-vector. */
} vd;
/// Vector of complex floating point numbers
typedef struct {
us size;
bool _foreign_data;
c* _data; /**< Pointer set if data storage is
intern. If this is set to zero, the
vector is a sub-vector. */
} vc;
/// Dense matrix of floating point values
typedef struct {
us n_rows;
@ -46,10 +30,20 @@ typedef struct {
c* _data;
} cmat;
#define setvecval(vec,index,val) \
dbgassert((((us) index) <= (vec)->size),OUTOFBOUNDSVEC); \
(vec)->_data[index] = val;
typedef dmat vd;
typedef cmat vc;
#define assert_equalsize(a,b) \
dbgassert((a)->n_rows == (b)->n_rows,SIZEINEQUAL); \
dbgassert((a)->n_cols == (b)->n_cols,SIZEINEQUAL);
#define is_vx(vx) ((vx)->n_cols == 1)
#define assert_vx(vx) dbgassert(is_vx(vx),"Not a vector!")
#define setvecval(vec,index,val) \
assert_vx(vec); \
dbgassert((((us) index) < (vec)->n_rows),OUTOFBOUNDSVEC); \
(vec)->_data[index] = val;
#define setmatval(mat,row,col,val) \
dbgassert((((us) row) <= mat->n_rows),OUTOFBOUNDSMATR); \
@ -62,22 +56,6 @@ typedef struct {
* @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 pointer to a value from a complex vector
*
* @param mat The vector
* @param row The row
*/
static inline c* getvcval(const vc* vec,us row){
dbgassert(row < vec->size,OUTOFBOUNDSVEC);
return &vec->_data[row];
}
/**
* Return a value from a matrix of floating points
*
@ -91,7 +69,6 @@ static inline d* getdmatval(const dmat* mat,us row,us col){
dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC);
return &mat->_data[col*mat->stride+row];
}
/**
* Return a value from a matrix of complex floating points
*
@ -106,17 +83,29 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){
return &mat->_data[col*mat->stride+row];
}
static inline d* getvdval(const vd* vec,us row){
dbgassert(vec,NULLPTRDEREF);
assert_vx(vec);
return getdmatval(vec,row,0);
}
/**
* Return pointer to a value from a complex vector
*
* @param mat The vector
* @param row The row
*/
static inline c* getvcval(const vc* vec,us row){
dbgassert(vec,NULLPTRDEREF);
assert_vx(vec);
return getcmatval(vec,row,0);
}
#ifdef LASP_DEBUG
#define OVERFLOW_MAGIC_NUMBER (-10e-45)
#define check_overflow_vx(vx) \
TRACE(15,"Checking overflow " #vx); \
if(!(vx)._foreign_data) { \
dbgassert((vx)._data[(vx).size] == OVERFLOW_MAGIC_NUMBER, \
"Buffer overflow detected on" #vx ); \
} \
#define check_overflow_xmat(xmat) \
TRACE(15,"Checking overflow " #xmat); \
if(!(xmat)._foreign_data) { \
@ -125,31 +114,13 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){
"Buffer overflow detected on" #xmat ); \
} \
#define check_overflow_vx check_overflow_xmat
#else
#define check_overflow_vx(vx)
#define check_overflow_xmat(xmat)
#endif
/**
* Sets all values in a vector to the value
*
* @param vec the vector to set
* @param value
*/
static inline void vd_set(vd* vec, d value){
d_set(vec->_data,value,vec->size);
}
/**
* 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){
c_set(vec->_data,value,vec->size);
}
/**
* Sets all values in a matrix to the value
*
@ -164,7 +135,7 @@ static inline void dmat_set(dmat* mat,const d value){
}
}
}
#define vd_set dmat_set
/**
* Sets all values in a matrix to the value
@ -180,49 +151,8 @@ static inline void cmat_set(cmat* mat,const c value){
}
}
}
#define vc_set cmat_set
/**
* Allocate data for a float vector.
*
* @param size Size of the vector
*
* @return vd with allocated data
*/
static inline vd vd_alloc(us size) {
vd result = { size, NULL,NULL};
#ifdef LASP_DEBUG
result._data = (d*) a_malloc((size+1)*sizeof(d));
result._data[size] = OVERFLOW_MAGIC_NUMBER;
#else
result._data = (d*) a_malloc(size*sizeof(d));
#endif // LASP_DEBUG
result._foreign_data = false;
#ifdef LASP_DEBUG
vd_set(&result,NAN);
#endif // LASP_DEBUG
return result;
}
/**
* Allocate data for a complex vector.
*
* @param size Size of the vector
*
* @return vc with allocated data
*/
static inline vc vc_alloc(us size) {
vc result = { size, NULL, NULL};
#ifdef LASP_DEBUG
result._data = (c*) a_malloc((size+1)*sizeof(c));
result._data[size] = OVERFLOW_MAGIC_NUMBER;
#else
result._data = (c*) a_malloc(size*sizeof(c));
#endif // LASP_DEBUG
result._foreign_data = false;
#ifdef LASP_DEBUG
vc_set(&result,NAN+I*NAN);
#endif // LASP_DEBUG
return result;
}
/**
* Allocate data for a matrix of floating points
*
@ -250,7 +180,6 @@ static inline dmat dmat_alloc(us n_rows,
return result;
}
/**
* Allocate a matrix of complex floating points
*
@ -276,6 +205,30 @@ static inline cmat cmat_alloc(const us n_rows,
#endif // LASP_DEBUG
return result;
}
/**
* Allocate data for a float vector.
*
* @param size Size of the vector
*
* @return vd with allocated data
*/
static inline vd vd_alloc(us size) {
return dmat_alloc(size,1);
}
/**
* Allocate data for a complex vector.
*
* @param size Size of the vector
*
* @return vc with allocated data
*/
static inline vc vc_alloc(us size) {
return cmat_alloc(size,1);
}
/**
* Creates a dmat from foreign data. Does not copy the data, but only
* initializes the row pointers. Assumes column-major ordering for the
@ -288,32 +241,62 @@ static inline cmat cmat_alloc(const us n_rows,
*
* @return
*/
static inline dmat dmat_foreign(const us n_rows,
const us n_cols,
d* data) {
dbgassert(data,NULLPTRDEREF);
dmat result = {n_rows,n_cols,true,n_rows,data};
return result;
}
static inline dmat dmat_foreign_vd(vd* vector) {
dbgassert(vector,NULLPTRDEREF);
dmat result = {vector->size,1,true,vector->size,vector->_data};
static inline dmat dmat_foreign(dmat* other) {
dbgassert(other,NULLPTRDEREF);
dmat result = {other->n_rows,
other->n_cols,
true,
other->stride,
other->_data};
return result;
}
/**
* Create vd from foreign data
* Create a dmat from foreign data. Assumes the stride of the data is
* n_rows.
*
* @param size Size of the vector
* @param data Pointer to data
* @param n_rows Number of rows
* @param n_cols Number of columns
* @param data Pointer to data storage
*
* @return
* @return dmat
*/
static inline vd vd_foreign(const us size,d* data) {
static inline dmat dmat_foreign_data(us n_rows,
us n_cols,
d* data,
bool own_data) {
dbgassert(data,NULLPTRDEREF);
vd result = {size,true,data};
dmat result = {n_rows,
n_cols,
!own_data,
n_rows,
data};
return result;
}
/**
* Create a cmat from foreign data. Assumes the stride of the data is
* n_rows.
*
* @param n_rows Number of rows
* @param n_cols Number of columns
* @param data Pointer to data storage
*
* @return dmat
*/
static inline cmat cmat_foreign_data(us n_rows,
us n_cols,
c* data,
bool own_data) {
dbgassert(data,NULLPTRDEREF);
cmat result = {n_rows,
n_cols,
!own_data,
n_rows,
data};
return result;
}
/**
* Creates a cmat from foreign data. Does not copy the data, but only
* initializes the row pointers. Assumes column-major ordering for the
@ -326,34 +309,18 @@ static inline vd vd_foreign(const us size,d* data) {
*
* @return
*/
static inline cmat cmat_foreign(const us n_rows,
const us n_cols,
c* data) {
dbgassert(data,NULLPTRDEREF);
cmat result = {n_rows,n_cols,true,n_rows,data};
static inline cmat cmat_foreign(cmat* other) {
dbgassert(other,NULLPTRDEREF);
cmat result = {other->n_rows,
other->n_cols,
true,
other->stride,
other->_data};
return result;
}
/**
* Free's data of a vector. Is safe to run on sub-vecs as well, to
* make API consistent. (Only free's data if data pointer is set)
*
* @param f Vector to free
*/
static inline void vd_free(vd* f) {
dbgassert(f,NULLPTRDEREF);
if(!(f->_foreign_data)) a_free(f->_data);
}
/**
* Free's data of a vector. Is safe to run on sub-vecs as well, to
* make API consistent. (Only free's data if data pointer is set)
*
* @param f Vector to free
*/
static inline void vc_free(vc* f) {
dbgassert(f,NULLPTRDEREF);
if(!(f->_foreign_data)) a_free(f->_data);
}
/**
* Free's data of dmat. Safe to run on sub-matrices as well.
*
@ -363,6 +330,8 @@ static inline void dmat_free(dmat* m) {
dbgassert(m,NULLPTRDEREF);
if(!(m->_foreign_data)) a_free(m->_data);
}
#define vd_free dmat_free
/**
* Free's data of dmat. Safe to run on sub-matrices as well.
*
@ -372,7 +341,7 @@ static inline void cmat_free(cmat* m) {
dbgassert(m,NULLPTRDEREF);
if(!(m->_foreign_data)) a_free(m->_data);
}
#define vc_free cmat_free
/**
* Copy some rows from one matrix to another
@ -457,40 +426,6 @@ static inline cmat cmat_submat(cmat* parent,
return result;
}
/**
* Copy contents of one vector to another
*
* @param to : Vector to write to
* @param from : Vector to read from
*/
static inline void vd_copy(vd* to,const vd* from) {
dbgassert(to && from,NULLPTRDEREF);
dbgassert(to->size==from->size,SIZEINEQUAL);
d_copy(to->_data,from->_data,to->size,1,1);
}
static inline void vd_copy_rows(vd* to,
const vd* from,
const us startrow_to,
const us startrow_from,
const us nrows) {
dbgassert(to && from,NULLPTRDEREF);
dbgassert(startrow_from+nrows <= from->size,OUTOFBOUNDSMATR);
dbgassert(startrow_to+nrows <= to->size,OUTOFBOUNDSMATR);
d_copy(&to->_data[startrow_to],
&from->_data[startrow_from],
nrows,1,1);
}
/**
* Copy contents of one vector to another
*
* @param to : Vector to write to
* @param from : Vector to read from
*/
static inline void vc_copy(vc* to,const vc* from) {
dbgassert(to && from,NULLPTRDEREF);
dbgassert(to->size==from->size,SIZEINEQUAL);
c_copy(to->_data,from->_data,to->size);
}
/**
* Copy contents of one matrix to another. Sizes should be equal
*
@ -524,6 +459,31 @@ static inline void cmat_copy(cmat* to,const cmat* from) {
}
}
/**
* Copy contents of one vector to another
*
* @param to : Vector to write to
* @param from : Vector to read from
*/
static inline void vd_copy(vd* to,const vd* from) {
dbgassert(to && from,NULLPTRDEREF);
assert_vx(to);
assert_vx(from);
dmat_copy(to,from);
}
/**
* Copy contents of one vector to another
*
* @param to : Vector to write to
* @param from : Vector to read from
*/
static inline void vc_copy(vc* to,const vc* from) {
dbgassert(to && from,NULLPTRDEREF);
assert_vx(to);
assert_vx(from);
cmat_copy(to,from);
}
/**
* Get a reference to a column of a matrix as a vector
@ -534,7 +494,7 @@ static inline void cmat_copy(cmat* to,const cmat* from) {
* @return vector with reference to column
*/
static inline vd dmat_column(dmat* x,us col) {
vd res = {x->n_rows,true,getdmatval(x,0,col)};
vd res = {x->n_rows,1,true,x->stride,getdmatval(x,0,col)};
return res;
}
@ -547,7 +507,7 @@ static inline vd dmat_column(dmat* x,us col) {
* @return vector with reference to column
*/
static inline vc cmat_column(cmat* x,us col) {
vc res = {x->n_rows,true,getcmatval(x,0,col)};
vc res = {x->n_rows,1,true,x->stride,getcmatval(x,0,col)};
return res;
}
@ -557,11 +517,14 @@ static inline vc cmat_column(cmat* x,us col) {
* @param a
* @param b
*/
static inline void vc_conj(vc* a,const vc* b) {
static inline void cmat_conj(cmat* a,const cmat* b) {
fsTRACE(15);
dbgassert(a && b,NULLPTRDEREF);
dbgassert(a->size == b->size,SIZEINEQUAL);
carray_conj(a->_data,b->_data,a->size);
dbgassert(a->n_cols == b->n_cols,SIZEINEQUAL);
dbgassert(a->n_rows == b->n_rows,SIZEINEQUAL);
for(us col=0;col<a->n_cols;col++) {
carray_conj(getcmatval(a,0,col),getcmatval(b,0,col),a->n_rows);
}
feTRACE(15);
}
@ -570,7 +533,7 @@ static inline void vc_conj(vc* a,const vc* b) {
*
* @param x
*/
static inline void cmat_conj(cmat* x) {
static inline void cmat_conj_inplace(cmat* x) {
dbgassert(x,NULLPTRDEREF);
for(us col=0;col<x->n_cols;col++) {
c_conj_inplace(getcmatval(x,0,col),x->n_rows);
@ -580,14 +543,15 @@ static inline void cmat_conj(cmat* x) {
#ifdef LASP_DEBUG
void print_cmat(const cmat* m);
void print_vc(const vc* m);
void print_vd(const vd* m);
void print_dmat(const dmat* m);
#define print_vc print_cmat
#define print_vd print_dmat
#else
#define print_cmat(m)
#define print_vc(m)
#define print_dmat(m)
#endif
#endif // LASP_MATH_H
#endif // LASP_MAT_H
//////////////////////////////////////////////////////////////////////

2
lasp/c/lasp_ps.c
View File

@ -148,7 +148,7 @@ void PowerSpectra_compute(const PowerSpectra* ps,
vc j_vec = cmat_column(&fft_work,j);
/* Compute the conjugate of spectra j */
vc_conj(&j_vec_conj,&j_vec);
cmat_conj(&j_vec_conj,&j_vec);
check_overflow_xmat(fft_work);

1
lasp/c/lasp_python.h
View File

@ -8,7 +8,6 @@
#pragma once
#ifndef LASP_PYTHON_H
#define LASP_PYTHON_H
#define TRACERPLUS (-5)
#include <numpy/ndarrayobject.h>
#ifdef LASP_DOUBLE_PRECISION
#define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT64

6
lasp/c/lasp_signals.h
View File

@ -8,7 +8,7 @@
#pragma once
#ifndef LASP_SIGNALS_H
#define LASP_SIGNALS_H
#include "lasp_math.h"
#include "lasp_mat.h"
/**
* Compute the signal power, that is \f$ \frac{1}{N} \sum_{i=0}^{N-1}
@ -19,10 +19,10 @@
*/
static inline d signal_power(vd* signal) {
d res = 0;
for(us i=0;i<signal->size;i++) {
for(us i=0;i<signal->n_rows;i++) {
res+= d_pow(*getvdval(signal,i),2);
}
res /= signal->size;
res /= signal->n_rows;
return res;
}

5
lasp/c/lasp_window.c
View File

@ -59,7 +59,10 @@ static d bartlett(us n,us N) {
}
int window_create(const WindowType wintype,vd* result,d* win_pow) {
fsTRACE(15);
us nfft = result->size;
dbgassert(result && win_pow,NULLPTRDEREF);
assert_vx(result);
us nfft = result->n_rows;
d (*win_fun)(us,us);
switch (wintype) {
case Hann: {

2
lasp/c/lasp_window.h
View File

@ -8,7 +8,7 @@
#pragma once
#ifndef LASP_WINDOW_H
#define LASP_WINDOW_H
#include "lasp_math.h"
#include "lasp_mat.h"
typedef enum {
Hann = 0,

23
lasp/config.pxi.in
View File

@ -29,7 +29,7 @@ cdef extern from "lasp_tracer.h":
void fsTRACE(int)
void feTRACE(int)
void clearScreen()
cdef extern from "lasp_math.h":
cdef extern from "lasp_mat.h":
ctypedef struct dmat:
us n_cols
us n_rows
@ -37,16 +37,17 @@ cdef extern from "lasp_math.h":
bint _foreign_data
ctypedef struct cmat:
pass
ctypedef struct vd:
pass
ctypedef struct vc:
pass
dmat dmat_foreign(us n_rows,
us n_cols,
d* data)
cmat cmat_foreign(us n_rows,
us n_cols,
c* data)
ctypedef cmat vc
ctypedef dmat vd
dmat dmat_foreign_data(us n_rows,
us n_cols,
d* data,
bint own_data)
cmat cmat_foreign_data(us n_rows,
us n_cols,
c* data,
bint own_data)
cmat cmat_alloc(us n_rows,us n_cols)
dmat dmat_alloc(us n_rows,us n_cols)

79
lasp/wrappers.pyx
View File

@ -50,13 +50,15 @@ cdef class Fft:
# result[:,:] = np.nan+1j*np.nan
cdef c[::1,:] result_view = result
cdef cmat r = cmat_foreign(result.shape[0],
result.shape[1],
&result_view[0,0])
cdef cmat r = cmat_foreign_data(result.shape[0],
result.shape[1],
&result_view[0,0],
False)
cdef dmat t = dmat_foreign(timedata.shape[0],
timedata.shape[1],
&timedata[0,0])
cdef dmat t = dmat_foreign_data(timedata.shape[0],
timedata.shape[1],
&timedata[0,0],
False)
Fft_fft(self._fft,&t,&r)
@ -74,18 +76,20 @@ cdef class Fft:
# result[:,:] = np.nan+1j*np.nan
cdef cmat f = cmat_foreign(freqdata.shape[0],
freqdata.shape[1],
&freqdata[0,0])
cdef cmat f = cmat_foreign_data(freqdata.shape[0],
freqdata.shape[1],
&freqdata[0,0],
False)
timedata = np.empty((nfft,nchannels),
dtype=NUMPY_FLOAT_TYPE,
order='F')
cdef d[::1,:] timedata_view = timedata
cdef dmat t = dmat_foreign(timedata.shape[0],
timedata.shape[1],
&timedata_view[0,0])
cdef dmat t = dmat_foreign_data(timedata.shape[0],
timedata.shape[1],
&timedata_view[0,0],
False)
Fft_ifft(self._fft,&f,&t)
@ -142,9 +146,10 @@ cdef class PowerSpectra:
cmat result_mat
td = dmat_foreign(nfft,
nchannels,
&timedata[0,0])
td = dmat_foreign_data(nfft,
nchannels,
&timedata[0,0],
False)
# The array here is created in such a way that the strides
# increase with increasing dimension. This is required for
@ -159,9 +164,10 @@ cdef class PowerSpectra:
cdef c[::1,:,:] result_view = result
result_mat = cmat_foreign(nfft//2+1,
nchannels*nchannels,
&result_view[0,0,0])
result_mat = cmat_foreign_data(nfft//2+1,
nchannels*nchannels,
&result_view[0,0,0],
False)
@ -208,8 +214,8 @@ cdef class AvPowerSpectra:
cdef vd weighting_vd
cdef vd* weighting_ptr = NULL
if(weighting.size != 0):
weighting_vd = vd_foreign(weighting.size,
&weighting[0])
weighting_vd = dmat_foreign_data(weighting.n_rows,1,
&weighting[0],False)
weighting_ptr = &weighting_vd
self.aps = AvPowerSpectra_alloc(nfft,
@ -242,9 +248,10 @@ cdef class AvPowerSpectra:
if nchannels != self.nchannels:
raise RuntimeError('Invalid number of channels')
td = dmat_foreign(nsamples,
nchannels,
&timedata[0,0])
td = dmat_foreign_data(nsamples,
nchannels,
&timedata[0,0],
False)
result_ptr = AvPowerSpectra_addTimeData(self.aps,
&td)
@ -260,9 +267,10 @@ cdef class AvPowerSpectra:
dtype = NUMPY_COMPLEX_TYPE,
order='F')
cdef c[::1,:,:] result_view = result
cdef cmat res = cmat_foreign(self.nfft//2+1,
nchannels*nchannels,
&result_view[0,0,0])
cdef cmat res = cmat_foreign_data(self.nfft//2+1,
nchannels*nchannels,
&result_view[0,0,0],
True)
# Copy result
cmat_copy(&res,result_ptr)
@ -282,7 +290,11 @@ cdef class FilterBank:
cdef:
c_FilterBank* fb
def __cinit__(self,d[::1,:] h, us nfft):
cdef dmat hmat = dmat_foreign(h.shape[0],h.shape[1],&h[0,0])
cdef dmat hmat = dmat_foreign_data(h.shape[0],
h.shape[1],
&h[0,0],
True)
self.fb = FilterBank_create(&hmat,nfft)
dmat_free(&hmat)
if not self.fb:
@ -293,7 +305,9 @@ cdef class FilterBank:
FilterBank_free(self.fb)
def filter_(self,d[:] input_):
cdef vd input_vd = vd_foreign(input_.size,&input_[0])
cdef dmat input_vd = dmat_foreign_data(input_.shape[0],1,
&input_[0],False)
cdef dmat output = FilterBank_filter(self.fb,&input_vd)
# Steal the pointer from output
@ -326,9 +340,10 @@ cdef class Decimator:
def decimate(self,d[::1,:] samples):
assert samples.shape[1] == self.nchannels,'Invalid number of channels'
cdef dmat d_samples = dmat_foreign(samples.shape[0],
samples.shape[1],
&samples[0,0])
cdef dmat d_samples = dmat_foreign_data(samples.shape[0],
samples.shape[1],
&samples[0,0],
False)
cdef dmat res = Decimator_decimate(self.dec,&d_samples)
result = dmat_to_ndarray(&res,True)
@ -337,4 +352,4 @@ cdef class Decimator:
def __dealloc__(self):
if self.dec != NULL:
Decimator_free(self.dec)
Decimator_free(self.dec)

2
test/test_bf.c
View File

@ -5,7 +5,7 @@
// Description:
//
//////////////////////////////////////////////////////////////////////
#include "lasp_math.h"
#include "lasp_mat.h"
int main() {

2
test/test_math.c
View File

@ -7,7 +7,7 @@
//////////////////////////////////////////////////////////////////////
#include "lasp_fft.h"
#include "lasp_math.h"
#include "lasp_mat.h"
#include "lasp_tracer.h"
#include "lasp_alg.h"