From f37ce6b5961d3bad33cc294fd7b3c249ea89bfb8 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong" Date: Sun, 4 Mar 2018 19:48:46 +0100 Subject: [PATCH] Fixed bug in filtebank, renamed ASCEE to LASP in lasp_python.h, added decimator. --- CMakeLists.txt | 3 +- lasp/c/CMakeLists.txt | 3 +- lasp/c/lasp_alg.c | 5 +- lasp/c/lasp_assert.h | 14 +- lasp/c/lasp_decimation.c | 214 ++++++++++++++++++++ lasp/c/lasp_decimation.h | 54 +++++ lasp/c/lasp_fft.c | 4 +- lasp/c/lasp_filterbank.c | 3 +- lasp/c/lasp_filterbank.h | 6 +- lasp/c/lasp_math.h | 8 +- lasp/c/lasp_math_raw.h | 10 +- lasp/c/lasp_python.h | 16 +- lasp/fir_design/octave_filter_design_fir.py | 27 +-- lasp/wrappers.pyx | 40 +++- test/filterbank_test.py | 2 +- 15 files changed, 359 insertions(+), 50 deletions(-) create mode 100644 lasp/c/lasp_decimation.c create mode 100644 lasp/c/lasp_decimation.h diff --git a/CMakeLists.txt b/CMakeLists.txt index ed76c55..834f2ee 100644 --- a/CMakeLists.txt +++ b/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" ) diff --git a/lasp/c/CMakeLists.txt b/lasp/c/CMakeLists.txt index 382b915..58cf81d 100644 --- a/lasp/c/CMakeLists.txt +++ b/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 ) diff --git a/lasp/c/lasp_alg.c b/lasp/c/lasp_alg.c index 2839792..1bd70a9 100644 --- a/lasp/c/lasp_alg.c +++ b/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;jn_cols;j++){ for(i=0;in_rows;i++) { - c* Aij = &A->data[i+j*n_rows]; - b->data[i] += *Aij * x->data[j]; + c* Aij = getcmatval(A,i,j); + b->_data[i] += *Aij * *getvcval(x,j); } diff --git a/lasp/c/lasp_assert.h b/lasp/c/lasp_assert.h index 82be4d1..e9584e9 100644 --- a/lasp/c/lasp_assert.h +++ b/lasp/c/lasp_assert.h @@ -6,8 +6,8 @@ // Basic tools for debugging using assert statements including text. ////////////////////////////////////////////////////////////////////// #pragma once -#ifndef ASCEE_ASSERT_H -#define ASCEE_ASSERT_H +#ifndef LASP_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 -#include "types.h" +#ifdef LASP_DEBUG +#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 ////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_decimation.c b/lasp/c/lasp_decimation.c new file mode 100644 index 0000000..af78ad6 --- /dev/null +++ b/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;channelnofbs[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;coln_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;chanfbs[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;channchannels;chan++) { + FilterBank_free(dec->fbs[chan]); + } + + a_free(dec->fbs); + a_free(dec); + + feTRACE(15); +} +////////////////////////////////////////////////////////////////////// + diff --git a/lasp/c/lasp_decimation.h b/lasp/c/lasp_decimation.h new file mode 100644 index 0000000..207a4dc --- /dev/null +++ b/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 +////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_fft.c b/lasp/c/lasp_fft.c index 5a956c5..162c0c9 100644 --- a/lasp/c/lasp_fft.c +++ b/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 */ diff --git a/lasp/c/lasp_filterbank.c b/lasp/c/lasp_filterbank.c index 92429cf..65b7b6c 100644 --- a/lasp/c/lasp_filterbank.c +++ b/lasp/c/lasp_filterbank.c @@ -127,6 +127,7 @@ dmat FilterBank_filter(FilterBank* fb, for(us col=0;colfilters,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); diff --git a/lasp/c/lasp_filterbank.h b/lasp/c/lasp_filterbank.h index 3d7a828..17f40bc 100644 --- a/lasp/c/lasp_filterbank.h +++ b/lasp/c/lasp_filterbank.h @@ -10,8 +10,8 @@ // each filter in the filterbank. ////////////////////////////////////////////////////////////////////// #pragma once -#ifndef FILTERBANK_H -#define FILTERBANK_H +#ifndef LASP_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 ////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_math.h b/lasp/c/lasp_math.h index f834d2b..dc03761 100644 --- a/lasp/c/lasp_math.h +++ b/lasp/c/lasp_math.h @@ -399,7 +399,7 @@ static inline void dmat_copy_rows(dmat* to,const dmat* from, for(col=0;colsize==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;coln_cols;col++) { d_copy(getdmatval(to,0,col), getdmatval(from,0,col), - to->n_rows); + to->n_rows,1,1); } } /** diff --git a/lasp/c/lasp_math_raw.h b/lasp/c/lasp_math_raw.h index 4a71708..3767b94 100644 --- a/lasp/c/lasp_math_raw.h +++ b/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 -#ifdef ASCEE_DOUBLE_PRECISION -#define ASCEE_NUMPY_FLOAT_TYPE NPY_FLOAT64 -#define ASCEE_NUMPY_COMPLEX_TYPE NPY_COMPLEX128 +#ifdef LASP_DOUBLE_PRECISION +#define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT64 +#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 ////////////////////////////////////////////////////////////////////// diff --git a/lasp/fir_design/octave_filter_design_fir.py b/lasp/fir_design/octave_filter_design_fir.py index 15c15a7..06da2d4 100644 --- a/lasp/fir_design/octave_filter_design_fir.py +++ b/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], --3:[4,8], --4:[4,8], --5:[4,8,4], --6:[4,8,4] +-2:[4,4], +-3:[4,4], +-4:[4,4,4], +-5:[4,4,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) diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index 81685ad..9e90049 100644 --- a/lasp/wrappers.pyx +++ b/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": - pass + ctypedef struct c_FilterBank "FilterBank" 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) \ No newline at end of file diff --git a/test/filterbank_test.py b/test/filterbank_test.py index 1d1c53b..79961df 100644 --- a/test/filterbank_test.py +++ b/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