Fixed bug in filtebank, renamed ASCEE to LASP in lasp_python.h, added decimator.

CMakeLists.txt

@@ -73,7 +73,8 @@ set(CYTHON_EXTRA_C_FLAGS "-Wno-sign-compare -Wno-cpp -Wno-implicit-fallthrough\
 
 # General make flags
 set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fPIC -std=c11  -Wall -Wextra -Wno-type-limits \
- -Werror=implicit-function-declaration -Werror=incompatible-pointer-types")
+ -Werror=implicit-function-declaration -Werror=incompatible-pointer-types \
+-Werror=return-type")
 
 # Debug make flags
 set(CMAKE_C_FLAGS_DEBUG "-g" )

lasp/c/CMakeLists.txt

@@ -16,8 +16,9 @@ add_library(lasp_lib
   lasp_mq.c
   lasp_worker.c
   lasp_dfifo.c
-  # lasp_filterbank.c
+  lasp_filterbank.c
   lasp_octave_fir.c
+  lasp_decimation.c
   )
 
 

lasp/c/lasp_alg.c

@@ -52,7 +52,6 @@ void cmv_dot(const cmat* A,const vc* restrict x,vc* restrict b){
 				
     #else
     size_t i,j;
-    size_t n_rows = A->n_rows;
 
     vc_set(b,0.0);
 
@@ -62,8 +61,8 @@ void cmv_dot(const cmat* A,const vc* restrict x,vc* restrict b){
     for(j=0;j<A->n_cols;j++){
         for(i=0;i<A->n_rows;i++) {
 
-            c* Aij = &A->data[i+j*n_rows];
+            c* Aij = getcmatval(A,i,j);
-            b->data[i] += *Aij * x->data[j];
+            b->_data[i] += *Aij * *getvcval(x,j);
 
         }
 

lasp/c/lasp_assert.h

@@ -6,8 +6,8 @@
 // Basic tools for debugging using assert statements including text.
 //////////////////////////////////////////////////////////////////////
 #pragma once
-#ifndef ASCEE_ASSERT_H
+#ifndef LASP_ASSERT_H
-#define ASCEE_ASSERT_H
+#define LASP_ASSERT_H
 
 #define OUTOFBOUNDSMATR "Out of bounds access on matrix row"
 #define OUTOFBOUNDSMATC "Out of bounds access on matrix column"
@@ -16,8 +16,8 @@
 #define ALLOCFAILED "Memory allocation failure in: "
 #define NULLPTRDEREF "Null pointer dereference in: "
 
-#ifdef ASCEE_DEBUG
+#ifdef LASP_DEBUG
-#include "types.h"
+#include "lasp_types.h"
 
 
 void DBG_AssertFailedExtImplementation(const char* file,
@@ -30,13 +30,13 @@ void DBG_AssertFailedExtImplementation(const char* file,
     	DBG_AssertFailedExtImplementation(__FILE__, __LINE__, assert_string );  \
     }
 
-#else // ASCEE_DEBUG not defined
+#else // LASP_DEBUG not defined
 
 #define dbgassert(assertion, assert_string)
 
 #define assertvalidptr(ptr) dbgassert(ptr,NULLPTRDEREF)
 
-#endif  // ASCEE_DEBUG
+#endif  // LASP_DEBUG
 
-#endif // ASCEE_ASSERT_H
+#endif // LASP_ASSERT_H
 //////////////////////////////////////////////////////////////////////

lasp/c/lasp_decimation.c

@@ -0,0 +1,214 @@
+// lasp_decimation.c
+//
+// Author: J.A. de Jong -ASCEE
+// 
+// Description:
+// Implementation of the decimator
+//////////////////////////////////////////////////////////////////////
+#include "lasp_decimation.h"
+#include "lasp_filterbank.h"
+#include "lasp_tracer.h"
+#include "lasp_alloc.h"
+#include "lasp_dfifo.h"
+
+#define DEC_FILTER_MAX_TAPS (128)
+#define DEC_FFT_LEN (1024)
+
+typedef struct {
+    DEC_FAC df;
+    us dec_fac;
+    us ntaps;
+    d h[DEC_FILTER_MAX_TAPS];
+} DecFilter;
+
+static __thread DecFilter DecFilters[] = {
+    {DEC_FAC_4,4,128,
+#include "dec_filter.c"     
+    }
+};
+
+
+typedef struct Decimator_s {
+    us nchannels;
+    us dec_fac;
+    FilterBank** fbs;
+    dFifo* output_fifo;
+} Decimator;
+
+
+Decimator* Decimator_create(us nchannels,DEC_FAC df) {
+    fsTRACE(15);
+
+    /* Find the filter */
+    const us nfilters = sizeof(DecFilters)/sizeof(DecFilter);
+    DecFilter* filter = DecFilters;
+    bool found = false;
+    for(us filterno = 0;filterno < nfilters; filterno++) {
+        if(filter->df == df) {
+            TRACE(15,"Found filter");
+            found = true;
+            break;
+        }
+        filter++;
+    }
+    if(!found) {
+        WARN("Decimation factor not found in list of filters");
+        return NULL;
+    }
+    
+    /* Create the filterbanks */
+    Decimator* dec = a_malloc(sizeof(Decimator));
+    dec->fbs = a_malloc(sizeof(FilterBank*)*nchannels);
+    dec->nchannels = nchannels;
+    dec->dec_fac = filter->dec_fac;
+
+    dmat h = dmat_foreign(filter->ntaps,1,filter->h);
+
+    for(us channelno=0;channelno<nchannels;channelno++) {
+        dec->fbs[channelno] = FilterBank_create(&h,DEC_FFT_LEN);
+    }
+
+    dmat_free(&h);
+    
+    /* Create input and output fifo's */
+    dec->output_fifo = dFifo_create(nchannels,DEC_FFT_LEN);
+
+    feTRACE(15);
+    return dec; 
+}
+
+static void lasp_downsample(dmat* downsampled,
+                            const dmat* filtered,
+                            us dec_fac) {
+    fsTRACE(15);
+    dbgassert(downsampled && filtered,NULLPTRDEREF);
+    dbgassert(filtered->n_rows/dec_fac == downsampled->n_rows,
+              "Incompatible number of rows");
+    dbgassert(downsampled->n_cols == filtered->n_cols,
+              "Incompatible number of rows");
+    
+    dbgassert(dec_fac> 0,"Invalid dec_fac");
+
+    for(us col=0;col<downsampled->n_cols;col++) {
+        /* Low-level BLAS copy. */
+        d_copy(getdmatval(downsampled,0,col),
+               getdmatval(filtered,0,col),
+               downsampled->n_rows,     /* number */
+               1,                       /* incx out */
+               dec_fac);                /* incx in */
+    }
+
+    check_overflow_xmat(*downsampled);
+    check_overflow_xmat(*filtered);
+
+    feTRACE(15);
+}
+
+dmat Decimator_decimate(Decimator* dec,const dmat* samples) {
+
+    fsTRACE(15);
+    dbgassert(dec && samples,NULLPTRDEREF);
+    const us nchannels = dec->nchannels;
+    const us dec_fac = dec->dec_fac;
+    dbgassert(samples->n_cols == nchannels,"Invalid number of "
+              "channels in samples");
+    
+
+    /* Not downsampled, but filtered result */
+    dmat filtered;
+
+    /* Filter each channel and store result in filtered. In first
+     * iteration the right size for filtered is allocated. */
+    for(us chan=0;chan<nchannels;chan++) {
+        vd samples_channel = dmat_column((dmat*)samples,
+                                         chan);
+
+        /* Low-pass filter stuff */
+        dmat filtered_res = FilterBank_filter(dec->fbs[chan],
+                                                 &samples_channel);
+
+        dbgassert(filtered_res.n_cols == 1,"Bug in FilterBank");
+
+        vd_free(&samples_channel);
+
+        if(chan==0) {
+            /* Allocate space for result */
+            filtered = dmat_alloc(filtered_res.n_rows,
+                                  nchannels);
+            
+        }
+        dmat filtered_col = dmat_submat(&filtered,
+                                        0,chan,
+                                        filtered_res.n_rows,
+                                        1);
+        
+        dbgassert(filtered_res.n_rows == filtered_col.n_rows,
+                  "Not all FilterBank's have same output number"
+                  " of rows!");
+        
+        dmat_copy_rows(&filtered_col,
+                       &filtered_res,
+                       0,0,filtered_res.n_rows);
+                                   
+        dmat_free(&filtered_res);
+        dmat_free(&filtered_col);
+        vd_free(&samples_channel);
+
+    }
+    
+    /* Push filtered result into output fifo */
+    dFifo_push(dec->output_fifo,
+               &filtered);
+
+    dmat_free(&filtered);
+
+
+    /* Now, downsample stuff */
+    dmat downsampled;
+    uVARTRACE(15,dec_fac);
+    us fifo_size = dFifo_size(dec->output_fifo);
+    if((fifo_size / dec_fac) > 0) {
+        
+        filtered = dmat_alloc((fifo_size/dec_fac)*dec_fac,
+                              nchannels);
+
+        int nsamples = dFifo_pop(dec->output_fifo,
+                                 &filtered,0);
+
+        dbgassert((us) nsamples % dec_fac == 0 && nsamples > 0,
+                  "BUG");
+        dbgassert((us) nsamples == (fifo_size/dec_fac)*dec_fac,"BUG");
+        
+        downsampled = dmat_alloc(nsamples/dec_fac,nchannels);
+        /* Do the downsampling work */
+        lasp_downsample(&downsampled,&filtered,dec_fac);
+        
+        dmat_free(&filtered);
+    }
+    else {
+        TRACE(15,"Empty downsampled");
+        downsampled = dmat_alloc(0,0);
+    }
+    
+    feTRACE(15);
+    /* return filtered; */
+    return downsampled;
+}
+
+
+void Decimator_free(Decimator* dec) {
+    fsTRACE(15);
+    dbgassert(dec,NULLPTRDEREF);
+    dFifo_free(dec->output_fifo);
+
+    for(us chan=0;chan<dec->nchannels;chan++) {
+        FilterBank_free(dec->fbs[chan]);
+    }
+
+    a_free(dec->fbs);
+    a_free(dec);
+
+    feTRACE(15);
+}
+//////////////////////////////////////////////////////////////////////
+

lasp/c/lasp_decimation.h

@@ -0,0 +1,54 @@
+// lasp_decimation.h
+//
+// Author: J.A. de Jong - ASCEE
+//
+// Description: Sample decimator, works on a whole block of
+// (uninterleaved) time samples at once. Decimates (downsamples and
+// low-pass filters by a factor given integer factor.
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef LASP_DECIMATION_H
+#define LASP_DECIMATION_H
+#include "lasp_types.h"
+#include "lasp_math.h"
+
+typedef struct Decimator_s Decimator;
+
+typedef enum DEC_FAC_e {
+    DEC_FAC_4 = 0,              // Decimate by a factor of 4
+} DEC_FAC;
+
+/** 
+ * Create a decimation filter for a given number of channels and
+ * decimation factor
+ *
+ * @param nchannels Number of channels
+
+ * @param d Decimation factor. Should be one of the implemented
+ * ones. (TODO: add more. At this point we have only a decimation
+ * factor of 4 implemented)
+ *
+ * @return Decimator handle. NULL on error
+ */
+Decimator* Decimator_create(us nchannels,DEC_FAC d);
+
+/** 
+ * Decimate given samples
+ *
+ * @param samples 
+ *
+ * @return Decimated samples. Can be an empty array.
+ */
+dmat Decimator_decimate(Decimator* dec,const dmat* samples);
+
+
+d Decimator_get_cutoff(Decimator*);
+/** 
+ * Free memory corresponding to Decimator
+ *
+ * @param dec Decimator handle.
+ */
+void Decimator_free(Decimator* dec);
+
+#endif // LASP_DECIMATION_H
+//////////////////////////////////////////////////////////////////////

lasp/c/lasp_fft.c

@@ -66,7 +66,7 @@ void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) {
     d* result_ptr = getvdval(result,0);
 
     /* Copy freqdata, to fft_result. */
-    d_copy(&result_ptr[1],&freqdata_ptr[2],nfft-1);
+    d_copy(&result_ptr[1],&freqdata_ptr[2],nfft-1,1,1);
     result_ptr[0] = freqdata_ptr[0];
 
     /* Perform inplace backward transform */
@@ -129,7 +129,7 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
      * part of the DC component equals zero. */
 
     /* Copy timedata, as it will be overwritten in the fft pass. */
-    d_copy(&result_ptr[1],getvdval(timedata,0),nfft);
+    d_copy(&result_ptr[1],getvdval(timedata,0),nfft,1,1);
 
 
     /* Perform fft */

lasp/c/lasp_filterbank.c

@@ -127,6 +127,7 @@ dmat FilterBank_filter(FilterBank* fb,
         for(us col=0;col<nfilters;col++) {
             vc output_fft_col = cmat_column(&output_fft,col);
             vc filter_col = cmat_column(&fb->filters,col);
+
             vc_hadamard(&output_fft_col,
                         &input_fft_col,
                         &filter_col);
@@ -139,7 +140,7 @@ dmat FilterBank_filter(FilterBank* fb,
         Fft_ifft(fb->fft,&output_fft,&output_block);
 
         dmat valid_samples = dmat_submat(&output_block,
-                                         0,0, /* startrow, startcol */
+                                         fb->P_m_1,0, /* startrow, startcol */
                                          nfft-fb->P_m_1, /* Number of rows */
                                          output_block.n_cols);
         dFifo_push(fb->output_fifo,&valid_samples);

lasp/c/lasp_filterbank.h

@@ -10,8 +10,8 @@
 // each filter in the filterbank.
 //////////////////////////////////////////////////////////////////////
 #pragma once
-#ifndef FILTERBANK_H
+#ifndef LASP_FILTERBANK_H
-#define FILTERBANK_H
+#define LASP_FILTERBANK_H
 #include "lasp_types.h"
 #include "lasp_math.h"
 typedef struct FilterBank_s FilterBank;
@@ -55,5 +55,5 @@ dmat FilterBank_filter(FilterBank* fb,
 void FilterBank_free(FilterBank* f);
 
 
-#endif // FILTERBANK_H
+#endif // LASP_FILTERBANK_H
 //////////////////////////////////////////////////////////////////////

lasp/c/lasp_math.h

@@ -399,7 +399,7 @@ static inline void dmat_copy_rows(dmat* to,const dmat* from,
     for(col=0;col<ncols;col++) {
         d* to_d = getdmatval(to,startrow_to,col);
         d* from_d = getdmatval(from,startrow_from,col);
-        d_copy(to_d,from_d,nrows);
+        d_copy(to_d,from_d,nrows,1,1);
     }
 }
 /** 
@@ -471,7 +471,7 @@ static inline cmat cmat_submat(cmat* parent,
 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);
+    d_copy(to->_data,from->_data,to->size,1,1);
 }
 static inline void vd_copy_rows(vd* to,
                                 const vd* from,
@@ -483,7 +483,7 @@ static inline void vd_copy_rows(vd* to,
     dbgassert(startrow_to+nrows <= to->size,OUTOFBOUNDSMATR);
     d_copy(&to->_data[startrow_to],
            &from->_data[startrow_from],
-           nrows);
+           nrows,1,1);
 }
 /** 
  * Copy contents of one vector to another
@@ -509,7 +509,7 @@ static inline void dmat_copy(dmat* to,const dmat* from) {
     for(us col=0;col<to->n_cols;col++) {
         d_copy(getdmatval(to,0,col),
                getdmatval(from,0,col),
-               to->n_rows);
+               to->n_rows,1,1);
     }
 }
 /** 

lasp/c/lasp_math_raw.h

@@ -174,13 +174,17 @@ static inline d d_dot(const d a[],const d b[],const us size){
  * @param from : Array to read from
  * @param size : Size of arrays
  */
-static inline void d_copy(d to[],const d from[],const us size){
+static inline void d_copy(d to[],
+                          const d from[],
+                          const us size,
+                          const us to_inc,
+                          const us from_inc){
     #if LASP_USE_BLAS == 1
-    cblas_dcopy(size,from,1,to,1);
+    cblas_dcopy(size,from,from_inc,to,to_inc);
     #else
     us i;
     for(i=0;i<size;i++)
-        to[i] = from[i];
+        to[i*to_inc] = from[i*from_inc];
     #endif
 }
 

lasp/c/lasp_python.h

@@ -6,14 +6,14 @@
 // Some routines to generate numpy arrays from matrices and vectors.
 //////////////////////////////////////////////////////////////////////
 #pragma once
-#ifndef ASCEE_PYTHON_H
+#ifndef LASP_PYTHON_H
-#define ASCEE_PYTHON_H
+#define LASP_PYTHON_H
 #include <numpy/ndarrayobject.h>
-#ifdef ASCEE_DOUBLE_PRECISION
+#ifdef LASP_DOUBLE_PRECISION
-#define ASCEE_NUMPY_FLOAT_TYPE NPY_FLOAT64
+#define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT64
-#define ASCEE_NUMPY_COMPLEX_TYPE NPY_COMPLEX128
+#define LASP_NUMPY_COMPLEX_TYPE NPY_COMPLEX128
 #else
-#define ASCEE_NUMPY_FLOAT_TYPE NPY_FLOAT32
+#define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT32
 #endif
 
 
@@ -28,7 +28,7 @@ static inline PyObject* dmat_to_ndarray(dmat* mat,bool transfer_ownership) {
     // more easy than using the PyArray_New syntax.
     npy_intp dims[] = {mat->n_cols,mat->n_rows};
     PyObject* arr_t = PyArray_SimpleNewFromData(2,dims,
-                                              ASCEE_NUMPY_FLOAT_TYPE,
+                                              LASP_NUMPY_FLOAT_TYPE,
                                               mat->_data);
     if(!arr_t) {
         WARN("Array creation failure");
@@ -56,5 +56,5 @@ static inline PyObject* dmat_to_ndarray(dmat* mat,bool transfer_ownership) {
 }
 
 
-#endif // ASCEE_PYTHON_H
+#endif // LASP_PYTHON_H
 //////////////////////////////////////////////////////////////////////

lasp/fir_design/octave_filter_design_fir.py

@@ -44,9 +44,9 @@ cut_fac_l = {
 0:0.98,
 -1:0.96,
 -2:.99,
--3:0.98,
+-3:0.96,
 -4:0.96,
--5:0.98,
+-5:0.96,
 -6:0.96
 }
 
@@ -59,9 +59,9 @@ cut_fac_u = {
 0:1.01,
 -1:1.02,
 -2:1.006,
--3:1.01,
+-3:1.02,
 -4:1.02,
--5:1.01,
+-5:1.02,
 -6:1.02
 }
 
@@ -73,11 +73,11 @@ decimation = {
 1:[4],
 0:[4],
 -1:[4],
--2:[4,8],
+-2:[4,4],
--3:[4,8],
+-3:[4,4],
--4:[4,8],
+-4:[4,4,4],
--5:[4,8,4],
+-5:[4,4,4],
--6:[4,8,4]
+-6:[4,4,4,4]
         }
 
 # Generate the header file
@@ -154,6 +154,8 @@ __thread OctaveFIR OctaveFIRs[%i] = {
         # Cut-off frequency of the filter
         f2 = fm*G**(1/(2*b))*cut_fac_u[x]
 
+        fc = fd/2/1.4
+
         fir = bandpass_fir_design(L,fd,f1,f2)        
         
         struct += "{ "
@@ -161,7 +163,7 @@ __thread OctaveFIR OctaveFIRs[%i] = {
             struct += "%0.16e," %i
         struct = struct[:-1] + " }},"
         
-        freq = np.logspace(np.log10(f1/3),np.log10(f2*3),1000)
+        freq = np.logspace(np.log10(f1/3),np.log10(fd),1000)
         H = freqResponse(fir,freq,fd)
         if showfig:        
             ax.semilogx(freq,20*np.log10(np.abs(H)))
@@ -173,8 +175,9 @@ __thread OctaveFIR OctaveFIRs[%i] = {
             ax.set_title('x = %i, fnom = %s' %(x,nominals[x]) )
         
             ax.set_ylim(-10,1)
-            ax.set_xlim(f1/1.1,f2*1.1)
+            ax.set_xlim(f1/1.1,fd)
-    
+            ax.axvline(fd/2)
+            ax.axvline(fc,color='red')
     struct+="\n};"
     cfile.write(struct)    
     

lasp/wrappers.pyx

@@ -226,8 +226,6 @@ cdef class AvPowerSpectra:
             dmat td
             cmat* result_ptr
 
-        if nsamples < self.nfft:
-            raise RuntimeError('Number of samples should be > nfft')
         if nchannels != self.nchannels:
             raise RuntimeError('Invalid number of channels')
 
@@ -261,8 +259,7 @@ cdef class AvPowerSpectra:
         return result
 
 cdef extern from "lasp_filterbank.h":
-    ctypedef struct c_FilterBank "FilterBank":
+    ctypedef struct c_FilterBank "FilterBank"
-        pass
     c_FilterBank* FilterBank_create(const dmat* h,const us nfft) nogil
     dmat FilterBank_filter(c_FilterBank* fb,const vd* x) nogil
     void FilterBank_free(c_FilterBank* fb) nogil
@@ -293,3 +290,38 @@ cdef class FilterBank:
         vd_free(&input_vd)
 
         return result
+
+cdef extern from "lasp_decimation.h":
+    ctypedef struct c_Decimator "Decimator"
+    ctypedef enum DEC_FAC:
+        DEC_FAC_4
+        
+    c_Decimator* Decimator_create(us nchannels,DEC_FAC d) nogil
+    dmat Decimator_decimate(c_Decimator* dec,const dmat* samples) nogil
+    void Decimator_free(c_Decimator* dec) nogil
+    
+cdef class Decimator:
+    cdef:
+        c_Decimator* dec
+        us nchannels
+    def __cinit__(self, us nchannels,us dec_fac):
+        assert dec_fac == 4, 'Invalid decimation factor'
+        self.nchannels = nchannels
+        self.dec = Decimator_create(nchannels,DEC_FAC_4)
+        if not self.dec:
+            raise RuntimeError('Error creating 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 res = Decimator_decimate(self.dec,&d_samples)
+        result = dmat_to_ndarray(&res,True)
+        dmat_free(&res)
+        return result
+        
+    def __dealloc__(self):
+        if self.dec != NULL:
+            Decimator_free(self.dec)

test/filterbank_test.py

@@ -6,7 +6,7 @@ Created on Mon Jan 15 19:45:33 2018
 @author: anne
 """
 import numpy as np
-from beamforming import FilterBank, cls
+from lasp import FilterBank, cls
 cls()
 
 nfft=50