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

@@ -1 +1,3 @@
 from .wrappers import *
+from .tools import *
+P_REF = 2e-5

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;

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;col<result->n_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.
  *

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
  *

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

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);
+
 /** 
  * 
  *

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");
 

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)

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)