Added dFifo queue for samples. Improved aps.c to handle all sizes of time data to be added

beamforming/c/CMakeLists.txt

@@ -15,6 +15,7 @@ add_library(beamforming_lib
   ps.c
   mq.c
   worker.c
+  dfifo.c
   )
 
 

beamforming/c/aps.c

@@ -9,13 +9,16 @@
 #include "aps.h"
 #include "ps.h"
 #include "ascee_alg.h"
+#include "dfifo.h"
+
+/* Multiplication factor for the maximum size of the fifo queue. This
+ * factor is multiplied by nfft to obtain the maximum size of the
+ * fifo. */
+#define FIFO_SIZE_MULT (5)
 
 typedef struct AvPowerSpectra_s {
-    us os;              /**< Counter set to the position where
-                         * the next time block should be taken
-                         * from */
     us nfft, nchannels;
-    us oo;
+    us overlap;                 /* Number of samples to overlap */
 
     us naverages;               /* Counter that counts the number of
                                  * averages taken for the computation
@@ -26,6 +29,8 @@ typedef struct AvPowerSpectra_s {
                                  * equal to nfft. */
 
 
+    dFifo* fifo;          /* Sample fifo storage */
+
     cmat ps_storage;           /**< Here we store the averaged
                                 * results for each Cross-power
                                 * spectra computed so far. */
@@ -39,6 +44,19 @@ typedef struct AvPowerSpectra_s {
     
 } AvPowerSpectra;
 
+void AvPowerSpectra_free(AvPowerSpectra* aps) {
+    fsTRACE(15);
+
+    PowerSpectra_free(aps->ps);
+    dFifo_free(aps->fifo);
+    dmat_free(&aps->buffer);
+    cmat_free(&aps->ps_storage);
+    cmat_free(&aps->ps_single);
+    a_free(aps);
+    
+    feTRACE(15);
+}
+
 AvPowerSpectra* AvPowerSpectra_alloc(const us nfft,
                                      const us nchannels,
                                      const d overlap_percentage,
@@ -47,28 +65,25 @@ AvPowerSpectra* AvPowerSpectra_alloc(const us nfft,
     fsTRACE(15);
     
     /* Check nfft */
-    if(nfft % 2 != 0 || nfft > ASCEE_MAX_NFFT) {
+    if(nfft==0 || nfft % 2 != 0 || nfft > ASCEE_MAX_NFFT) {
-        WARN("nfft should be even");
+        WARN("Invalid nfft");
         feTRACE(15);
         return NULL;
     }
 
     /* Check overlap percentage */
-    if(overlap_percentage >= 100) {
+    if(overlap_percentage < 0. || overlap_percentage >= 100.) {
-        WARN("Overlap percentage >= 100!");
+        WARN("Invalid overlap percentage");
-        feTRACE(15);
-        return NULL;
-    }
-    if(overlap_percentage < 0) {
-        WARN("Overlap percentage should be positive!");
         feTRACE(15);
         return NULL;
     }
     
     /* Compute and check overlap offset */
-    us oo = (us) (((d) nfft)-overlap_percentage*((d) nfft)/100);
+    us overlap = (us) (overlap_percentage*((d) nfft)/100);
-    iVARTRACE(15,oo);
+    if(overlap == nfft) {
-    if(oo == 0) {oo++;}
+        WARN("Overlap percentage results in full overlap, decreasing overlap.");
+        overlap--;
+    }
     
     PowerSpectra* ps = PowerSpectra_alloc(nfft,wt);
     if(!ps) {
@@ -83,13 +98,14 @@ AvPowerSpectra* AvPowerSpectra_alloc(const us nfft,
     aps->nfft = nfft;
     aps->ps = ps;
     aps->naverages = 0;
-    aps->oo = oo;
+    aps->overlap = overlap;
-    aps->os = oo;
+
+    aps->buffer = dmat_alloc(nfft,nchannels);
 
     /* Allocate vectors and matrices */
-    aps->buffer = dmat_alloc(nfft,nchannels);
     aps->ps_storage = cmat_alloc(nfft/2+1,nchannels*nchannels);
     aps->ps_single = cmat_alloc(nfft/2+1,nchannels*nchannels);
+    aps->fifo = dFifo_create(nchannels,FIFO_SIZE_MULT*nfft);
     
     cmat_set(&aps->ps_storage,0);
     feTRACE(15);
@@ -151,92 +167,65 @@ cmat* AvPowerSpectra_addTimeData(AvPowerSpectra* aps,
     const us nfft = aps->nfft;
     
     dbgassert(timedata->n_cols == nchannels,"Invalid time data");
-    dbgassert(timedata->n_rows >= nfft,"Invalid time data. "
+    dbgassert(timedata->n_rows > 0,"Invalid time data. "
-              "Should at least have nfft rows");
+              "Should at least have one row");
 
-    const us oo = aps->oo;
+    const us nsamples = timedata->n_rows;
-    us* os = &aps->os;
 
-    us os_timedata = 0;
+    /* Split up timedata in blocks of size ~ (FIFO_SIZE_MULT-1)nfft */
+    const us max_blocksize = (FIFO_SIZE_MULT-1)*nfft;
 
-    dmat buffer = aps->buffer;
+    us pos = 0; /* Current position in timedata buffer */
 
-    /* Retrieve the buffer and use it to make the first time block. */
+    /* dFifo handle */
-    if(*os < oo) {
+    dFifo* fifo = aps->fifo;
-        TRACE(15,"Using saved data from previous run");
-        dbgassert(false,"not tested")
-        dmat tmp = dmat_alloc(nfft,nchannels);
-        dbgassert(0 <= *os,"BUG");
-        dbgassert(*os <= nfft,"BUG"); 
 
-        /* copy_dmat_rows(&tmp, */
+    do {
-        /*                &buffer, */
+        us nsamples_part = pos+max_blocksize <= nsamples ?
-        /*                *os,       /\* Startrow_from *\/ */
+            max_blocksize : nsamples-pos;
-        /*                0,         /\* Startrow to *\/ */
-        /*                nfft - *os /\* nrows *\/ */
-        /*     ); */
 
-        /* copy_dmat_rows(&tmp, */
+        /* Obtain sub matrix */
-        /*                timedata, */
+        dmat timedata_part = dmat_submat(timedata,
-        /*                0, */
+                                         pos, /* Startrow */
-        /*                nfft - *os, */
+                                         0,   /* Startcol */
-        /*                *os */
+                                         nsamples_part, /* n_rows */
-        /*     ); */
+                                         nchannels);    /* n_cols */
     
-
+        if(dFifo_push(fifo,&timedata_part)!=SUCCESS) {
-        AvPowerSpectra_addBlock(aps,&tmp);
+            WARN("Fifo push failed.");
-
-        os_timedata = oo + *os - nfft;
-        dbgassert(os_timedata < nfft,"BUG");
-        dmat_free(&tmp);
         }
 
-    /* Run until we cannot go any further */
+        /* Temporary storage buffer */
-    while ((os_timedata + nfft) <= timedata->n_rows) {
+        dmat* buffer = &aps->buffer;
 
-        dmat tmp = dmat_submat(timedata,
+        /* Pop samples from the fifo while there are still at
-                               os_timedata, /* Startrow */
+         * least nfft samples available */
-                               0,           /* Start column */
+        while (dFifo_size(fifo) >= nfft) {
-                               nfft,        /* Number of rows */
+            int popped = dFifo_pop(fifo,
-                               nchannels);  /* Number of columns */
+                                   buffer,
+                                   aps->overlap); /* Keep 'overlap'
+                                                   * number of samples
+                                                   * in the queue */
 
+            dbgassert((us) popped == nfft,"Bug in dFifo");
             /* Process the block of time data */
-        AvPowerSpectra_addBlock(aps,&tmp);
+            AvPowerSpectra_addBlock(aps,buffer);
 
-        iVARTRACE(15,os_timedata);
-        os_timedata += oo;
-
-        dmat_free(&tmp);
-        iVARTRACE(15,os_timedata);
         }
 
-    /* We copy the last piece of samples from the timedata to the
+        dmat_free(&timedata_part);
-     * buffer */
-    dmat_copy_rows(&buffer,
-                   timedata,
-                   0,           /* startrow_to */
-                   timedata->n_rows-nfft, /* startrow_from */
-                   nfft);       /* Number of rows */
 
-    *os = os_timedata+nfft-timedata->n_rows;
+        /* Update position */
+        pos+=nsamples_part;
         
-    dbgassert(*os <= nfft,"BUG");
+    } while (pos < nsamples);
+    dbgassert(pos == nsamples,"BUG");
     
     feTRACE(15);
     return &aps->ps_storage;
 }
 
-void AvPowerSpectra_free(AvPowerSpectra* aps) {
-    fsTRACE(15);
 
-    PowerSpectra_free(aps->ps);
-    dmat_free(&aps->buffer);
-    cmat_free(&aps->ps_storage);
-    cmat_free(&aps->ps_single);
-    a_free(aps);
-    
-    feTRACE(15);
-}
 
 
 //////////////////////////////////////////////////////////////////////

beamforming/c/ascee_math.h

@@ -116,7 +116,7 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){
                   "Buffer overflow detected on" #vx );          \
     }                                                           \
     else {                                                      \
-    WARN("Cannot check overflow on foreign buffer");            \
+    DBGWARN("Cannot check overflow on foreign buffer");            \
     }
 
 #define check_overflow_xmat(xmat)                               \
@@ -127,7 +127,7 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){
                   "Buffer overflow detected on" #xmat );                \
     }                                                                   \
     else {                                                      \
-    WARN("Cannot check overflow on foreign buffer");            \
+    DBGWARN("Cannot check overflow on foreign buffer");            \
     }
 
 #else
@@ -401,7 +401,10 @@ static inline dmat dmat_submat(const dmat* parent,
     dbgassert(n_rows+startrow <= parent->n_rows,OUTOFBOUNDSMATR);
     dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC);
 
-    dmat result = { n_rows,n_cols,true,n_rows-startrow,
+    dmat result = { n_rows,n_cols,
+                    true,           // Foreign data = true
+                    parent->n_rows, // This is the stride to get to
+                    // the next column.
                     getdmatval(parent,startrow,startcol)};
 
     return result;
@@ -427,7 +430,11 @@ static inline cmat cmat_submat(cmat* parent,
     dbgassert(n_rows+startrow <= parent->n_rows,OUTOFBOUNDSMATR);
     dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC);
 
-    cmat result = { n_rows,n_cols,true,n_rows-startrow,
+
+    cmat result = { n_rows,n_cols,
+                    true,           // Foreign data = true
+                    parent->n_rows, // This is the stride to get to
+                    // the next column.
                     getcmatval(parent,startrow,startcol)};
 
     return result;

beamforming/c/ascee_math_raw.c

@@ -5,6 +5,7 @@
 // Description:
 // Operations working on raw arrays of floating point numbers
 //////////////////////////////////////////////////////////////////////
+#define TRACERPLUS (-5)
 #include "ascee_math_raw.h"
 #if ASCEE_USE_BLAS
 #include <cblas.h>

beamforming/c/dfifo.c

@@ -0,0 +1,138 @@
+// dfifo.c
+//
+// Author: J.A. de Jong -ASCEE
+// 
+// Description:
+// Implementation of the dFifo queue
+//////////////////////////////////////////////////////////////////////
+
+#include "dfifo.h"
+#define DFIFO_QUEUE_MAX_BLOCKS (50)
+
+typedef struct dFifo_s {
+    dmat queue;
+    us start_row;
+    us end_row;
+} dFifo;
+us dFifo_size(dFifo* fifo) {
+    fsTRACE(15);
+    dbgassert(fifo,NULLPTRDEREF);
+    dbgassert(fifo->start_row <= fifo->end_row,"BUG");
+    feTRACE(15);
+    
+    return fifo->end_row-fifo->start_row;
+}
+dFifo* dFifo_create(const us nchannels,
+                    const us max_size) {
+
+    fsTRACE(15);
+    dFifo* fifo = a_malloc(sizeof(dFifo));
+    fifo->queue = dmat_alloc(max_size,nchannels);
+    fifo->start_row = 0;
+    fifo->end_row = 0;
+    feTRACE(15);
+    return fifo;
+}
+void dFifo_free(dFifo* fifo) {
+    fsTRACE(15);
+    dmat_free(&fifo->queue);
+    a_free(fifo);
+    feTRACE(15);
+}
+int dFifo_push(dFifo* fifo,const dmat* data) {
+    fsTRACE(15);
+    dbgassert(fifo && data, NULLPTRDEREF);
+    dbgassert(data->n_cols == fifo->queue.n_cols,
+              "Invalid number of columns in data");
+
+
+    dmat queue = fifo->queue;
+    const us max_size = queue.n_rows;
+    const us nchannels = queue.n_cols;
+
+    us* start_row = &fifo->start_row;
+    us* end_row = &fifo->end_row;
+    const us added_size = data->n_rows;
+    const us size_before = dFifo_size(fifo);
+
+    if(added_size + dFifo_size(fifo) > max_size) {
+        return FIFO_QUEUE_FULL;
+    }
+
+    if(*end_row + added_size > max_size) {
+        if(size_before != 0) {
+            /* Shift the samples to the front of the queue. TODO: this
+             * might not be the most optimal implementation (but it is the
+             * most simple). */
+            TRACE(15,"Shift samples back in buffer");
+            uVARTRACE(15,size_before);
+            dmat tmp = dmat_alloc(size_before,nchannels);
+            TRACE(15,"SFSG");
+            dmat_copy_rows(&tmp,
+                           &queue,
+                           0,
+                           *start_row,
+                           size_before);
+            TRACE(15,"SFSG");
+            dmat_copy_rows(&queue,
+                           &tmp,
+                           0,0,size_before);
+
+            *end_row -= *start_row;
+            *start_row = 0;
+
+            dmat_free(&tmp);
+        }
+        else {
+            *start_row = 0;
+            *end_row = 0;
+        }
+    }
+
+    /* Now, copy samples */
+    dmat_copy_rows(&queue,      /* to */
+                   data,        /* from */
+                   *end_row, /* startrow_to */
+                   0,                /* startrow_from */
+                   added_size);      /* n_rows */
+
+    /* Increase the size */
+    *end_row += added_size;
+
+    feTRACE(15);
+    return SUCCESS;
+}
+int dFifo_pop(dFifo* fifo,dmat* data,const us keep) {
+    fsTRACE(15);
+    dbgassert(fifo && data,NULLPTRDEREF);
+    dbgassert(data->n_cols == fifo->queue.n_cols,
+              "Invalid number of columns in data");
+    dbgassert(keep < data->n_rows, "Number of samples to keep should"
+              " be smaller than requested number of samples");
+
+    us* start_row = &fifo->start_row;
+    us* end_row = &fifo->end_row;
+    us cur_contents = dFifo_size(fifo);
+    us requested = data->n_rows;
+    
+    us obtained = requested > cur_contents ? cur_contents : requested;
+    dbgassert(obtained > keep,"Number of samples to keep should be"
+              " smaller than requested number of samples");
+
+    uVARTRACE(15,requested);
+    uVARTRACE(15,obtained);
+    uVARTRACE(15,*start_row);
+    uVARTRACE(15,*end_row);
+    
+    dmat_copy_rows(data,
+                   &fifo->queue,
+                   0,
+                   *start_row,
+                   obtained);
+    
+    *start_row += obtained - keep;
+
+    feTRACE(15);
+    return (int) obtained;
+}
+//////////////////////////////////////////////////////////////////////

beamforming/c/dfifo.h

@@ -0,0 +1,76 @@
+// dfifo.h
+//
+// Author: J.A. de Jong - ASCEE
+//
+// Description:
+// API of a contiguous fifo buffer of samples.
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef DFIFO_H
+#define DFIFO_H
+#include "types.h"
+#include "ascee_math.h"
+typedef struct dFifo_s dFifo;
+
+/** 
+ * Create a fifo buffer
+ *
+ * @param nchannels Number of channels to store for
+ * @param max_size Maximum size of the queue.
+ *
+ * @return Pointer to fifo queue.
+ */
+dFifo* dFifo_create(const us nchannels,
+                    const us max_size);
+
+
+#define FIFO_QUEUE_FULL (-1)
+
+/** 
+ * Pushes samples into the fifo.
+ *
+ * @param fifo dFifo handle
+ *
+ * @param data data to push. Number of columns should be equal to
+ * nchannels.
+ *
+ * @return 0 on success, FIFO_QUEUE_FULL when samples do not fit.
+ */
+int dFifo_push(dFifo* fifo,const dmat* data);
+
+/** 
+ * Pop samples from the queue
+ *
+ * @param[in] fifo dFifo handle
+
+ * @param[out] data Pointer to dmat where popped data will be
+ * stored. Should have nchannels number of columns. If n_rows is
+ * larger than current storage, the queue is emptied.
+
+ * @param[in] keep Keeps a number of samples for the next dFifo_pop(). If
+ * keep=0, then no samples will be left. Keep should be smaller than
+ * the number of rows in data.
+ *
+ * @return Number of samples obtained in data.
+ */
+int dFifo_pop(dFifo* fifo,dmat* data,const us keep);
+
+/** 
+ * Returns current size of the fifo
+ *
+ * @param[in] fifo dFifo handle
+ *
+ * @return Current size
+ */
+us dFifo_size(dFifo* fifo);
+
+/** 
+ * Free a dFifo object
+ *
+ * @param[in] fifo dFifo handle.
+ */
+void dFifo_free(dFifo* fifo);
+
+#endif // DFIFO_H
+//////////////////////////////////////////////////////////////////////
+

beamforming/c/ps.c

@@ -144,12 +144,9 @@ void PowerSpectra_compute(const PowerSpectra* ps,
              * is documented */
             vc res = cmat_column(result,i+j*nchannels);
 
-            check_overflow_vx(res);
             vc i_vec = cmat_column(&fft_work,i);
             vc j_vec = cmat_column(&fft_work,j);
 
-            check_overflow_xmat(fft_work);
-
             /* Compute the conjugate of spectra j */
             vc_conj(&j_vec_conj,&j_vec);
 

beamforming/wrappers.pyx

@@ -1,6 +1,6 @@
 include "config.pxi"
 
-setTracerLevel(-5)
+setTracerLevel(15)
 cdef extern from "cblas.h":
     int openblas_get_num_threads()
     void openblas_set_num_threads(int)