First work on single pole lowpass IIR filter

lasp/c/CMakeLists.txt

@@ -19,7 +19,9 @@ add_library(lasp_lib
   lasp_filterbank.c
   lasp_octave_fir.c
   lasp_decimation.c
+  lasp_sp_lowpass.c
   )
 
 
+
 target_link_libraries(lasp_lib fftpack openblas)

lasp/c/lasp_sp_lowpass.c

@@ -0,0 +1,89 @@
+// lasp_sp_lowpass.c
+//
+// Author: J.A. de Jong -ASCEE
+// 
+// Description:
+// Single pole lowpass IIR filter implementation.
+//////////////////////////////////////////////////////////////////////
+#include "lasp_sp_lowpass.h"
+#include "lasp_alloc.h"
+
+typedef struct SPLowpass_s {
+
+    d a;
+    d b;
+
+    d xlast,ylast;
+
+} SPLowpass;
+
+SPLowpass* SPLowpass_create(d fs,d tau) {
+    fsTRACE(15);
+
+    d tau_i = 1/tau;
+    d T = 1/fs;
+
+    if(fs <= 0) {
+        WARN("Invalid sampling frequency");
+        return NULL;
+    }
+    else if(T >= tau ) {
+        WARN("Invalid tau, should be (much) larger than sampling"
+             " time (1/fs)");
+        return NULL;
+    }
+
+
+    SPLowpass* lp = a_malloc(sizeof(SPLowpass));
+    lp->xlast = 0;
+    lp->ylast = 0;
+
+    lp->a = tau_i/(2*fs+tau_i);
+    lp->b = (2*fs-tau_i)/(2*fs+tau_i);
+
+    feTRACE(15);
+    return lp;
+}
+
+vd SPLowpass_filter(SPLowpass* lp,
+                    const vd* input) {
+
+    fsTRACE(15);
+    dbgassert(lp && input,NULLPTRDEREF);
+    us input_size = input->size;
+    if(input_size == 0) {
+        return vd_alloc(0);
+    }
+    
+    vd output = vd_alloc(input_size);
+    const d xlast = lp->xlast;
+    const d ylast = lp->ylast;
+    const d a = lp->a;
+    const d b = lp->b;
+    
+    *getvdval(&output,0) = a*(xlast + *getvdval(input,0)) +
+        b*ylast;
+    
+    for(us i=1;i<input_size;i++) {
+
+        *getvdval(&output,i) = a*(*getvdval(input,i-1) +
+                                  *getvdval(input,i))  +
+            b*(*getvdval(&output,i-1));
+
+    }
+
+    lp->xlast = *getvdval(input  ,input_size-1);
+    lp->ylast = *getvdval(&output,input_size-1);
+
+    feTRACE(15);
+    return output;
+}
+
+void SPLowpass_free(SPLowpass* lp) {
+    fsTRACE(15);
+    dbgassert(lp,NULLPTRDEREF);
+    a_free(lp);
+    feTRACE(15);
+}
+
+//////////////////////////////////////////////////////////////////////

lasp/c/lasp_sp_lowpass.h

@@ -0,0 +1,51 @@
+// lasp_sp_lowpass.h
+//
+// Author: J.A. de Jong - ASCEE
+//
+// Description: Single-pole low pass IIR filter. Used for example for
+// fast and fast and slow averaging of the square of the squared
+// A-weighted acoustic pressure. This class uses the bilinear
+// transform to derive the filter coefficients.
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef LASP_SP_LOWPASS_H
+#define LASP_SP_LOWPASS_H
+#include "lasp_types.h"
+#include "lasp_math.h"
+
+typedef struct SPLowpass_s SPLowpass;
+
+/** 
+ * Create a single pole lowpass IIR filter.
+ *
+ * @param fs Sampling frequency [Hz]
+
+ * @param tau Time constant of the lowpass filter. Should be >0,
+ * otherwise an error occurs
+ *
+ * @return Pointer to dynamically allocated lowpass filter. NULL on
+ * error.
+ */
+SPLowpass* SPLowpass_create(d fs,d tau);
+
+/** 
+ * Use the filter to filter input data
+ *
+ * @param lp Lowpass filter handlen
+ * @param input Vector of input samples.
+ *
+ * @return Output data. Length is equal to input at all cases.
+ */
+vd SPLowpass_filter(SPLowpass* lp,
+                    const vd* input);
+
+/** 
+ * Free a single pole lowpass filter
+ *
+ * @param lp 
+ */
+void SPLowpass_free(SPLowpass* lp);
+
+
+#endif // LASP_SP_LOWPASS_H
+//////////////////////////////////////////////////////////////////////

lasp/wrappers.pyx

@@ -337,4 +337,37 @@ cdef class Decimator:
         
     def __dealloc__(self):
         if self.dec != NULL:
-            Decimator_free(self.dec)
+            Decimator_free(self.dec)
+
+
+cdef extern from "lasp_sp_lowpass.h":
+    ctypedef struct c_SPLowpass "SPLowpass"
+    c_SPLowpass* SPLowpass_create(d fs,d tau)
+    vd SPLowpass_filter(c_SPLowpass* lp,
+                        const vd* _input)
+    void SPLowpass_free(c_SPLowpass* lp)
+
+cdef class SPLowpass:
+    cdef:
+        c_SPLowpass* lp
+    def __cinit__(self,d fs,d tau):
+        self.lp = SPLowpass_create(fs,tau)
+        if not self.lp:
+            raise RuntimeError('Error creating lowpass filter')
+
+    def __dealloc__(self):
+        if self.lp:
+            SPLowpass_free(self.lp)
+
+    def filter_(self,d[:] input_):
+        # cdef vd input_vd = vd_foreign(input_.size,&input_[0])
+        # cdef dmat output = FilterBank_filter(self.fb,&input_vd)
+
+        # # Steal the pointer from output
+        # result = dmat_to_ndarray(&output,True)
+
+        # dmat_free(&output)
+        # vd_free(&input_vd)
+
+        # return result
+        pass