Bugfix in equalizer. Seems to be working now. Bugfix in wrapper to make equalizer return array one-dimensional

lasp/c/CMakeLists.txt

@@ -1,6 +1,7 @@
 if(!LASP_DEBUG)
 SET_SOURCE_FILES_PROPERTIES(lasp_sosfilterbank.c PROPERTIES COMPILE_FLAGS -O3)
 SET_SOURCE_FILES_PROPERTIES(lasp_slm.c PROPERTIES COMPILE_FLAGS -O3)
+SET_SOURCE_FILES_PROPERTIES(lasp_eq.c PROPERTIES COMPILE_FLAGS -O3)
 endif(!LASP_DEBUG)
 
 add_library(lasp_lib
@@ -21,6 +22,7 @@ add_library(lasp_lib
   lasp_sosfilterbank.c
   lasp_decimation.c
   lasp_slm.c
+  lasp_eq.c
   )
 
 

lasp/c/lasp_eq.c

@@ -0,0 +1,70 @@
+#include "lasp_eq.h"
+#include "lasp_assert.h"
+
+typedef struct Eq {
+    Sosfilterbank* fb;
+    us nfilters;
+    vd ampl_values;
+} Eq;
+
+Eq* Eq_create(Sosfilterbank* fb) {
+    fsTRACE(15);
+    assertvalidptr(fb);
+    Eq* eq = a_malloc(sizeof(Eq));
+    eq->fb = fb;
+    eq->nfilters = Sosfilterbank_getFilterbankSize(fb);
+    eq->ampl_values = vd_alloc(eq->nfilters);
+    vd_set(&(eq->ampl_values), 1.0);
+    feTRACE(15);
+    return eq;
+}
+vd Eq_equalize(Eq* eq,const vd* input_data) {
+    fsTRACE(15);
+    assertvalidptr(eq);
+    assert_vx(input_data);
+    vd result = vd_alloc(input_data->n_rows);
+    dmat_set(&result, 0);
+    dmat filtered = Sosfilterbank_filter(eq->fb, input_data);
+
+    for(us filter=0;filter<eq->nfilters;filter++) {
+        d ampl = *getvdval(&(eq->ampl_values), filter);
+        /// TODO: Replace this code with something more fast from BLAS.
+        for(us sample=0;sample<filtered.n_rows;sample++) {
+            d* res = getvdval(&result, sample);
+            *res = *res + *getdmatval(&filtered, sample, filter) * ampl;
+        }
+    }
+    dmat_free(&filtered);
+
+    feTRACE(15);
+    return result;
+}
+
+void Eq_setLevels(Eq* eq,const vd* levels) {
+    fsTRACE(15);
+    assertvalidptr(eq);
+    assert_vx(levels);
+    dbgassert(levels->n_rows == eq->nfilters, "Invalid levels size");
+    for(us ch=0;ch<eq->nfilters;ch++){
+        d level = *getvdval(levels, ch);
+        *getvdval(&(eq->ampl_values), ch) = d_pow(10, level/20);
+    }
+
+    feTRACE(15);
+}
+
+us Eq_getNLevels(const Eq* eq) {
+    fsTRACE(15);
+    assertvalidptr(eq);
+    feTRACE(15);
+    return eq->nfilters;
+}
+
+void Eq_free(Eq* eq) {
+    fsTRACE(15);
+    assertvalidptr(eq);
+    assertvalidptr(eq->fb);
+    Sosfilterbank_free(eq->fb);
+    vd_free(&(eq->ampl_values));
+    feTRACE(15);
+}

lasp/c/lasp_eq.h

@@ -0,0 +1,60 @@
+// lasp_eq.h
+//
+// Author: J.A. de Jong - ASCEE
+//
+// Description: Implementation of an equalizer using the Second Order Sections
+// filter bank implementation. 
+#pragma once
+#ifndef LASP_EQ_H
+#define LASP_EQ_H
+#include "lasp_sosfilterbank.h"
+typedef struct Eq Eq;
+
+/**
+ * Initialize an equalizer using the given Filterbank. Note: takes over pointer
+ * ownership of fb! Sets levels of all filterbanks initially to 0 dB.
+ *
+ * @param[in] fb Filterbank to be used in equalizer
+ * @return Equalizer handle, NULL on error
+ * */
+Eq* Eq_create(Sosfilterbank* fb);
+
+/**
+ * Equalize a given piece of data using current settings of the levels.
+ *
+ * @param[in] eq Equalizer handle
+ * @param[in] input_data Input data to equalize
+ * @return Equalized data. Newly allocated vector. Ownership is transferred.
+ * */
+vd Eq_equalize(Eq* eq,const vd* input_data);
+
+/**
+ * Returns number of channels of the equalizer. Note: takes over pointer
+ * ownership of fb!
+ *
+ * @param[in] eq Equalizer handle
+ * */
+us Eq_getNLevels(const Eq* eq);
+
+/**
+ * Set amplification values for each filter in the equalizer.
+ *
+ * @param[in] eq Equalizer handle
+ * @param[in] levels: Vector with level values for each channel. Should have
+ * length equal to the number of filters in the filterbank.
+ * */
+void Eq_setLevels(Eq* eq, const vd* levels);
+
+/** 
+ * Cleans up an existing Equalizer
+ *
+ * @param[in] eq Equalizer handle
+ */
+void Eq_free(Eq* eq);
+
+    
+        
+
+#endif // LASP_EQ_H
+// //////////////////////////////////////////////////////////////////////
+

lasp/c/lasp_siggen.c

@@ -5,7 +5,7 @@
 // Description:
 // Signal generator implementation
 //////////////////////////////////////////////////////////////////////
-#define TRACERPLUS (-5)
+/* #define TRACERPLUS (-5) */
 #include "lasp_siggen.h"
 #include "lasp_alloc.h"
 #include "lasp_assert.h"

lasp/filter/filterbank_design.py

@@ -466,8 +466,14 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
     def nominal_txt(self, x):
         # Put the nominal frequencies in a dictionary for easy access with
         # x as the key.
-        index = x - self.xs[0]
-        return self._nominal_txt[index]
+        if type(x) == int:
+            index = x - self.xs[0]
+            return self._nominal_txt[index]
+        elif type(x) == list:
+            index_start = x[0] - self.xs[0]
+            index_stop = x[-1] - self.xs[0]
+            return self._nominal_txt[index_start:index_stop+1]
+            
 
     def band_limits(self, x, filter_class=0):
         """Returns the third octave band filter limits for filter designator x.

lasp/lasp_octavefilter.py

@@ -228,6 +228,11 @@ class SosFilterBank:
             xmin: Minimum value for the bands
             xmax: Maximum value for the bands
         """
+        if xmin is None:
+            xmin = self.designer.xs[0]
+        if xmax is None:
+            xmax = self.designer.xs[-1]
+
         self.fs = fs
         self.xs = list(range(xmin, xmax + 1))
         nfilt = len(self.xs)
@@ -290,10 +295,6 @@ class SosThirdOctaveFilterBank(SosFilterBank):
             xmax: Maximum value for the bands
         """
         self.designer = ThirdOctaveBankDesigner(fs)
-        if xmin is None:
-            xmin = self.designer.xs[0]
-        if xmax is None:
-            xmax = self.designer.xs[-1]
         SosFilterBank.__init__(self, fs, xmin, xmax)
 
 
@@ -313,10 +314,6 @@ class SosOctaveFilterBank(SosFilterBank):
             xmax: Maximum value for the bands, if not specified, use maximum
         """
         self.designer = OctaveBankDesigner(fs)
-        if xmin is None:
-            xmin = self.designer.xs[0]
-        if xmax is None:
-            xmax = self.designer.xs[-1]
         SosFilterBank.__init__(self, fs, xmin, xmax)
 
 

lasp/lasp_slm.py

@@ -149,7 +149,7 @@ class SLM:
                                        'data': levels[:, [i]],
                                        'x': x}
         if self.include_overall and self.fbdesigner is not None:
-            output['overall'] = {'t': t, 'data': levels[:, i+1], 'x': 0}
+            output['overall'] = {'t': t, 'data': levels[:, [i+1]], 'x': 0}
 
         return output
 

lasp/wrappers.pyx

@@ -674,3 +674,58 @@ cdef class Siggen:
         siggen = Siggen()
         siggen._siggen = c_siggen
         return siggen
+
+cdef extern from "lasp_eq.h" nogil:
+    ctypedef struct c_Eq "Eq"
+    c_Eq* Eq_create(c_Sosfilterbank* fb)
+    vd Eq_equalize(c_Eq* eq,const vd* input_data) 
+    us Eq_getNLevels(const c_Eq* eq)
+    void Eq_setLevels(c_Eq* eq, const vd* levels)
+    void Eq_free(c_Eq* eq)
+
+cdef class Equalizer:
+    cdef:
+        c_Eq* ceq
+
+    def __cinit__(self, SosFilterBank cdef_fb):
+        """
+        Initialize equalizer using given filterbank. Note: Steals pointer of
+        underlying c_Sosfilterbank!!
+        """
+        self.ceq = Eq_create(cdef_fb.fb)
+        # Set this pointer to NULL, such that the underlying c_SosfilterBank is
+        # not deallocated.
+        cdef_fb.fb = NULL
+
+    def getNLevels(self):
+        return Eq_getNLevels(self.ceq)
+
+    def setLevels(self,d[:] new_levels):
+        cdef dmat dmat_new_levels = dmat_foreign_data(new_levels.shape[0],
+                                        1,
+                                        &new_levels[0],
+                                        False)
+        Eq_setLevels(self.ceq, &dmat_new_levels) 
+        dmat_free(&dmat_new_levels)
+
+    def equalize(self, d[::1] input_data):
+        cdef:
+            vd res
+        cdef dmat input_dmat = dmat_foreign_data(input_data.size,
+                                        1,
+                                        &input_data[0],
+                                        False)
+        with nogil:
+            res = Eq_equalize(self.ceq, &input_dmat)
+
+        # Steal the pointer from output
+        py_res = dmat_to_ndarray(&res,True)[:,0]
+
+        dmat_free(&res)
+        vd_free(&input_dmat)
+
+        return py_res
+
+    def __dealloc__(self):
+        if self.ceq:
+            Eq_free(self.ceq)