From b31866974e97e613b2e7b3f45c73d3bdd4b30021 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong" Date: Mon, 12 Feb 2018 19:32:23 +0100 Subject: [PATCH] Renamed copy_dmat_rows to dmat_copy_rows and changed API. Added P_REF. Added Hadamard product of two matrices. Added vd_copy_rows. Put some dynamic allocation of fft.c to object data. Bugfix in dbgassert ps.c. Moved getFreq to a tools.py --- beamforming/__init__.py | 2 ++ beamforming/c/aps.c | 6 ++---- beamforming/c/ascee_alg.h | 30 ++++++++++++++++++++++++++++++ beamforming/c/ascee_math.h | 16 ++++++++++++++-- beamforming/c/fft.c | 12 ++++++++---- beamforming/c/fft.h | 9 +++++++++ beamforming/c/ps.c | 5 +---- beamforming/tools.py | 13 +++++++++++++ beamforming/wrappers.pyx | 7 ++----- 9 files changed, 81 insertions(+), 19 deletions(-) create mode 100644 beamforming/tools.py diff --git a/beamforming/__init__.py b/beamforming/__init__.py index fff7b7d..e96bd63 100644 --- a/beamforming/__init__.py +++ b/beamforming/__init__.py @@ -1 +1,3 @@ from .wrappers import * +from .tools import * +P_REF = 2e-5 diff --git a/beamforming/c/aps.c b/beamforming/c/aps.c index f1c579c..375cf7b 100644 --- a/beamforming/c/aps.c +++ b/beamforming/c/aps.c @@ -128,13 +128,11 @@ static void AvPowerSpectra_addBlock(AvPowerSpectra* aps, cmat* ps_storage = &aps->ps_storage; c naverages = (++aps->naverages); - cVARTRACE(15,naverages); /* Scale previous result */ cmat_scale(ps_storage, (naverages-1)/naverages); - uVARTRACE(15,(us) aps->ps); PowerSpectra_compute(aps->ps, block, @@ -220,10 +218,10 @@ cmat* AvPowerSpectra_addTimeData(AvPowerSpectra* aps, /* We copy the last piece of samples from the timedata to the * buffer */ - copy_dmat_rows(&buffer, + dmat_copy_rows(&buffer, timedata, - timedata->n_rows-nfft, /* startrow_from */ 0, /* startrow_to */ + timedata->n_rows-nfft, /* startrow_from */ nfft); /* Number of rows */ *os = os_timedata+nfft-timedata->n_rows; diff --git a/beamforming/c/ascee_alg.h b/beamforming/c/ascee_alg.h index 9cfd54d..439b800 100644 --- a/beamforming/c/ascee_alg.h +++ b/beamforming/c/ascee_alg.h @@ -133,6 +133,36 @@ static inline void vc_hadamard(vc* result,const vc* a,const vc* b) { check_overflow_vx(*b); feTRACE(15); } +/** + * Compute the element-wise (Hadamard) product of a and b, and store + * in result + * + * @param[out] result + * @param[in] a + * @param[in] b + */ +static inline void cmat_hadamard(cmat* result, + const cmat* a, + const cmat* b) { + fsTRACE(15); + dbgassert(result && a && b,NULLPTRDEREF); + + dbgassert(result->n_rows==a->n_rows,SIZEINEQUAL); + dbgassert(result->n_cols==a->n_cols,SIZEINEQUAL); + dbgassert(b->n_rows==a->n_rows,SIZEINEQUAL); + dbgassert(b->n_cols==a->n_cols,SIZEINEQUAL); + + for(us col=0;coln_cols;col++) { + c_hadamard(result->col_ptrs[col], + a->col_ptrs[col], + b->col_ptrs[col], + a->n_rows); + } + + feTRACE(15); +} + + /** * Compute the matrix vector product for complex-valued types: b = A*x. * diff --git a/beamforming/c/ascee_math.h b/beamforming/c/ascee_math.h index 45838b0..ce524d9 100644 --- a/beamforming/c/ascee_math.h +++ b/beamforming/c/ascee_math.h @@ -397,9 +397,9 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){ * @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, +static inline void dmat_copy_rows(dmat* to,const dmat* from, us startrow_to, + us startrow_from, us nrows) { us col,ncols = to->n_cols; @@ -488,6 +488,18 @@ static inline void vd_copy(vd* to,const vd* from) { dbgassert(to->size==from->size,SIZEINEQUAL); d_copy(to->ptr,from->ptr,to->size); } +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->ptr[startrow_to], + &from->ptr[startrow_from], + nrows); +} /** * Copy contents of one vector to another * diff --git a/beamforming/c/fft.c b/beamforming/c/fft.c index 249a894..f39179e 100644 --- a/beamforming/c/fft.c +++ b/beamforming/c/fft.c @@ -14,6 +14,8 @@ typedef struct Fft_s { us nfft; vd fft_work; + vd fft_result; /**< Temporary storage for the FFT + * result */ } Fft; Fft* Fft_alloc(const us nfft) { @@ -29,6 +31,8 @@ Fft* Fft_alloc(const us nfft) { /* Initialize foreign fft lib */ fft->fft_work = vd_alloc(2*nfft+15); + fft->fft_result = vd_alloc(nfft); + npy_rffti(nfft,fft->fft_work.ptr); check_overflow_vx(fft->fft_work); @@ -41,6 +45,7 @@ void Fft_free(Fft* fft) { fsTRACE(15); dbgassert(fft,NULLPTRDEREF); vd_free(&fft->fft_work); + vd_free(&fft->fft_result); a_free(fft); feTRACE(15); } @@ -51,7 +56,6 @@ 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, "Invalid size for time data rows." " Should be equal to nfft"); @@ -59,8 +63,10 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) { dbgassert(result->size == (nfft/2+1),"Invalid number of rows in" " result array"); - vd fft_result = vd_alloc(nfft); + /* Obtain fft_result */ + vd fft_result = fft->fft_result; + /* Copy timedata, as it will be overwritten in the fft pass. */ vd_copy(&fft_result,timedata); /* Perform fft */ @@ -83,7 +89,6 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) { check_overflow_vx(fft_result); check_overflow_vx(fft->fft_work); - vd_free(&fft_result); feTRACE(15); } @@ -110,7 +115,6 @@ void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result) { check_overflow_xmat(*timedata); check_overflow_xmat(*result); - feTRACE(15); } diff --git a/beamforming/c/fft.h b/beamforming/c/fft.h index 5a8f98b..33bcf1e 100644 --- a/beamforming/c/fft.h +++ b/beamforming/c/fft.h @@ -57,6 +57,15 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result); */ void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result); +/** + * Perform inverse FFT + * + * @param[in] fft Fft handle + * @param[in] freqdata Frequency domain data + * @param[out] timedata + */ +void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata); + /** * * diff --git a/beamforming/c/ps.c b/beamforming/c/ps.c index 7e613ea..21c5d4c 100644 --- a/beamforming/c/ps.c +++ b/beamforming/c/ps.c @@ -76,16 +76,13 @@ void PowerSpectra_compute(const PowerSpectra* ps, dbgassert(ps && timedata && result,NULLPTRDEREF); - const us nchannels = timedata->n_rows; + const us nchannels = timedata->n_cols; const us nfft = Fft_nfft(ps->fft); uVARTRACE(15,nchannels); const d win_pow = ps->win_pow; dVARTRACE(15,win_pow); /* Sanity checks for the matrices */ - dbgassert(timedata->n_cols == nchannels,"timedata n_cols " - "should be equal to nchannels"); - dbgassert(timedata->n_rows == nfft,"timedata n_rows " "should be equal to nfft"); diff --git a/beamforming/tools.py b/beamforming/tools.py new file mode 100644 index 0000000..fbe8729 --- /dev/null +++ b/beamforming/tools.py @@ -0,0 +1,13 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""! +Author: J.A. de Jong - ASCEE + +Description: +""" +import numpy as np + +def getFreq(fs, nfft): + df = fs/nfft # frequency resolution + K = nfft//2+1 # number of frequency bins + return np.linspace(0, (K-1)*df, K) diff --git a/beamforming/wrappers.pyx b/beamforming/wrappers.pyx index 85ce3de..5250cc4 100644 --- a/beamforming/wrappers.pyx +++ b/beamforming/wrappers.pyx @@ -22,6 +22,8 @@ cdef extern from "fft.h": void Fft_fft(c_Fft*,dmat * timedate,cmat * res) nogil us Fft_nfft(c_Fft*) + + cdef class Fft: cdef: c_Fft* _fft @@ -227,8 +229,3 @@ cdef class AvPowerSpectra: dmat_free(&td) return result - - def getFreq(self, d fs): - cdef d df = fs/self.nfft # frequency resolution - cdef us K = self.nfft//2+1 # number of frequency bins - return np.linspace(0, (K-1)*df, K)