Partially done implementation of Octave and one-third octave filterbank.

lasp/c/dec_filter_4.c

@@ -0,0 +1,128 @@
+{0.000000000000000000e+00,
+    -2.185738981122740943e-06,
+    -2.794005403671438005e-06,
+    9.474117030816240549e-06,
+    4.047683164318276243e-05,
+    8.159721257084205745e-05,
+    1.081272986403995632e-04,
+    8.819754747279720866e-05,
+    1.292089375486244570e-18,
+    -1.501119890529942897e-04,
+    -3.184064926464186956e-04,
+    -4.310893684346957895e-04,
+    -4.094689334832080671e-04,
+    -2.059289266182427725e-04,
+    1.632175892397718904e-04,
+    6.039349831684049600e-04,
+    9.609848668279323634e-04,
+    1.065386827404308026e-03,
+    8.003236325878445483e-04,
+    1.624425171981582980e-04,
+    -7.076130570048387320e-04,
+    -1.544654803296158230e-03,
+    -2.031572983126828796e-03,
+    -1.908002710205690703e-03,
+    -1.080539502548179872e-03,
+    3.071493710103571440e-04,
+    1.878468916273352439e-03,
+    3.114730098982527052e-03,
+    3.515061533711236544e-03,
+    2.779296151548805299e-03,
+    9.515559570652546090e-04,
+    -1.533803816637046430e-03,
+    -3.936604506597823384e-03,
+    -5.417051906137869584e-03,
+    -5.310758368959681876e-03,
+    -3.387308238052255758e-03,
+    -2.020189585846376504e-17,
+    3.941575036075541973e-03,
+    7.193441582498504364e-03,
+    8.546931560120834062e-03,
+    7.242589074154118407e-03,
+    3.295234447391589671e-03,
+    -2.391580800508445910e-03,
+    -8.190127435210622919e-03,
+    -1.217754814426674249e-02,
+    -1.272615168819673029e-02,
+    -9.085393051846636994e-03,
+    -1.766149937181874050e-03,
+    7.423829658638344230e-03,
+    1.575380997299343638e-02,
+    2.029368152419728719e-02,
+    1.881247426920701696e-02,
+    1.060259142796622991e-02,
+    -3.026237731635000611e-03,
+    -1.877030360461057201e-02,
+    -3.192876928283294724e-02,
+    -3.747119479720732033e-02,
+    -3.132999667548151679e-02,
+    -1.158792905301173765e-02,
+    2.076675004965774715e-02,
+    6.175096553249778686e-02,
+    1.050334155991766993e-01,
+    1.431963133103141828e-01,
+    1.693251638559906402e-01,
+    1.785441122547395953e-01,
+    1.691180120289761668e-01,
+    1.428459416149487626e-01,
+    1.046477658307891495e-01,
+    6.144841153891738433e-02,
+    2.063940827275027520e-02,
+    -1.150252035146304662e-02,
+    -3.106003822342199780e-02,
+    -3.710127967927729503e-02,
+    -3.157313813741360886e-02,
+    -1.853722959527762462e-02,
+    -2.984746593388511431e-03,
+    1.044334023900628412e-02,
+    1.850493210242899755e-02,
+    1.993456916585750749e-02,
+    1.545344532020422393e-02,
+    7.271929013423104535e-03,
+    -1.727500581384247688e-03,
+    -8.873383043962101632e-03,
+    -1.241029376178266405e-02,
+    -1.185679908045850044e-02,
+    -7.961640663049178793e-03,
+    -2.321032949510159950e-03,
+    3.192604628017779635e-03,
+    7.004731963719125487e-03,
+    8.251271964876608425e-03,
+    6.931579083478471744e-03,
+    3.790698415906206074e-03,
+    -1.938928732007791559e-17,
+    -3.244201578800122790e-03,
+    -5.075193378062196545e-03,
+    -5.164847959311490676e-03,
+    -3.744259045513907962e-03,
+    -1.455155744093486990e-03,
+    9.003467998605640694e-04,
+    2.622284972952067337e-03,
+    3.306540383383136175e-03,
+    2.920618605100662978e-03,
+    1.755414747876200702e-03,
+    2.859847735621343567e-04,
+    -1.002153516667671955e-03,
+    -1.762137666026810951e-03,
+    -1.867706509456178859e-03,
+    -1.413021131661008033e-03,
+    -6.438025163091647607e-04,
+    1.469135513812169238e-04,
+    7.190488353391707912e-04,
+    9.501796827231699877e-04,
+    8.500266304384044257e-04,
+    5.292428644636570966e-04,
+    1.415165912868260117e-04,
+    -1.763689026852747436e-04,
+    -3.456944177857647653e-04,
+    -3.578084944368887165e-04,
+    -2.589144427740892114e-04,
+    -1.190202164875697825e-04,
+    9.922942950529402484e-19,
+    6.497054701652780085e-05,
+    7.525727688871557033e-05,
+    5.231825144807937484e-05,
+    2.280076622100770170e-05,
+    4.215016359691300044e-06,
+    -6.989289499831618075e-07,
+    -0.000000000000000000e+00}

lasp/c/lasp_decimation.c

@@ -11,7 +11,10 @@
 #include "lasp_alloc.h"
 #include "lasp_dfifo.h"
 
+// The Maximum number of taps in a decimation filter
 #define DEC_FILTER_MAX_TAPS (128)
+
+// The FFT length
 #define DEC_FFT_LEN (1024)
 
 typedef struct {
@@ -23,7 +26,7 @@ typedef struct {
 
 static __thread DecFilter DecFilters[] = {
     {DEC_FAC_4,4,128,
-#include "dec_filter.c"     
+#include "dec_filter_4.c"
     }
 };
 
@@ -112,13 +115,14 @@ dmat Decimator_decimate(Decimator* dec,const dmat* samples) {
     const us dec_fac = dec->dec_fac;
     dbgassert(samples->n_cols == nchannels,"Invalid number of "
               "channels in samples");
-    
+    dbgassert(samples->n_rows > 0,"Number of rows should be >0")
 
     /* 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);
@@ -131,35 +135,38 @@ dmat Decimator_decimate(Decimator* dec,const dmat* samples) {
 
         vd_free(&samples_channel);
 
-        if(chan==0) {
+        if(filtered_res.n_rows > 0) {
-            /* Allocate space for result */
+            if(chan==0) {
-            filtered = dmat_alloc(filtered_res.n_rows,
+                /* Allocate space for result */
-                                  nchannels);
+                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);
+        }
+        else {
+            filtered = dmat_alloc(0, 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);
-
     }
-    
+    if(filtered.n_rows > 0) {
-    /* Push filtered result into output fifo */
+        /* Push filtered result into output fifo */
-    dFifo_push(dec->output_fifo,
+        dFifo_push(dec->output_fifo,
-               &filtered);
+                   &filtered);
-
+    }
     dmat_free(&filtered);
 
 

lasp/c/lasp_decimation.h

@@ -15,7 +15,7 @@
 typedef struct Decimator_s Decimator;
 
 typedef enum DEC_FAC_e {
-    DEC_FAC_4 = 0,              // Decimate by a factor of 4
+  DEC_FAC_4 = 0, // Decimate by a factor of 4
 } DEC_FAC;
 
 /**

lasp/c/lasp_filterbank.c

@@ -56,6 +56,12 @@ FilterBank* FilterBank_create(const dmat* h,
     fb->output_fifo = dFifo_create(nfilters,FIFO_SIZE_MULT*nfft);
     fb->input_fifo = dFifo_create(1,FIFO_SIZE_MULT*nfft);
 
+    // Initialize the input fifo with zeros.
+    // dmat init_zero = dmat_alloc(nfft - P, 1);
+    // dmat_set(&init_zero,0);
+    // dFifo_push(fb->input_fifo, &init_zero);
+    // dmat_free(&init_zero);
+
     /* Create a temporary buffer which is going to be FFT'th to
      * contain the filter transfer functions.
      */

lasp/c/lasp_mat.h

@@ -129,7 +129,7 @@ static inline c* getvcval(const vc* vec,us row){
  */
 static inline void dmat_set(dmat* mat,const d value){
     dbgassert(mat,NULLPTRDEREF);
-    if(likely(mat->n_cols && mat->n_rows)) {
+    if(likely(mat->n_cols * mat->n_rows > 0)) {
         for(us col=0;col<mat->n_cols;col++) {
             d_set(getdmatval(mat,0,col),value,mat->n_rows);
         }
@@ -145,7 +145,7 @@ static inline void dmat_set(dmat* mat,const d value){
  */
 static inline void cmat_set(cmat* mat,const c value){
     dbgassert(mat,NULLPTRDEREF);
-    if(likely(mat->n_cols && mat->n_rows)) {
+    if(likely(mat->n_cols * mat->n_rows > 0)) {
         for(us col=0;col<mat->n_cols;col++) {
             c_set(getcmatval(mat,0,col),value,mat->n_rows);
         }

lasp/filter/bandpass_fir.py

@@ -0,0 +1,221 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE
+
+Description: FIR filter design for octave bands from 16Hz to 16 kHz for a
+sampling frequency of 48 kHz, FIR filter design for one-third octave bands
+
+See test/octave_fir_test.py for a testing
+
+"""
+from .fir_design import bandpass_fir_design, freqResponse as frsp
+import numpy as np
+__all__ = ['OctaveBankDesigner', 'ThirdOctaveBankDesigner']
+
+
+class FilterBankDesigner:
+    """
+    A class responsible for designing FIR filters
+    """
+    G = 10**(3/10)
+    fr = 1000.
+    L = 256  # Filter order
+    fs = 48000.  # Sampling frequency
+
+    def fm(self, x):
+        """
+        Returns the exact midband frequency of the bandpass filter
+
+        Args:
+            x: Midband designator
+        """
+        # Exact midband frequency
+        return self.G**(x/self.b)*self.fr
+
+    def fl(self, x):
+        """
+        Returns the exact cut-on frequency of the bandpass filter
+
+        Args:
+            x: Midband designator
+        """
+        return self.fm(x)*self.G**(-1/(2*self.b))
+
+    def fu(self, x):
+        """
+        Returns the exact cut-off frequency of the bandpass filter
+
+        Args:
+            x: Midband designator
+        """
+        return self.fm(x)*self.G**(1/(2*self.b))
+
+    def createFilter(self, fs, x):
+        """
+        Create a FIR filter for band designator b and sampling frequency fs.
+        Decimation should be obtained from decimation() method.
+
+        Returns:
+            filter: 1D ndarray with FIR filter coefficients
+        """
+        assert np.isclose(fs, self.fs), "Invalid sampling frequency"
+        fd = fs / np.prod(self.decimation(x))
+
+        # For designing the filter, the lower and upper frequencies need to be
+        # slightly adjusted to fall within the limits for a class 1 filter.
+        fl = self.fl(x)*self.fac_l(x)
+        fu = self.fu(x)*self.fac_u(x)
+
+        return bandpass_fir_design(self.L, fd, fl, fu)
+
+    def freqResponse(self, fs, x, freq):
+        """
+        Compute the frequency response for a certain filter
+
+        Args:
+            fs: Sampling frequency [Hz]
+            x: Midband designator
+        """
+        fir = self.createFilter(fs, x)
+
+        # Decimated sampling frequency [Hz]
+        fd = fs / np.prod(self.decimation(x))
+
+        return frsp(fd, freq, fir)
+
+
+class OctaveBankDesigner(FilterBankDesigner):
+    """
+    Octave band filter designer
+    """
+    def __init__(self):
+        pass
+
+    @property
+    def b(self):
+        # Band division, 1 for octave bands
+        return 1
+
+    @property
+    def xs(self):
+        return list(range(-6, 5))
+
+    def nominal(self, x):
+        # Text corresponding to the nominal frequency
+        nominals = {4: '16k',
+                    3: '8k',
+                    2: '4k',
+                    1: '2k',
+                    0: '1k',
+                    -1: '500',
+                    -2: '250',
+                    -3: '125',
+                    -4: '63',
+                    -5: '31.5',
+                    -6: '16'}
+        return nominals[x]
+
+    def fac_l(self, x):
+        """
+        Factor with which to multiply the cut-on frequency of the FIR filter
+        """
+        if x == 4:
+            return .995
+        elif x in (3, 1):
+            return .99
+        elif x in(-6, -4, -2, 2, 0):
+            return .98
+        else:
+            return .96
+
+    def fac_u(self, x):
+        """
+        Factor with which to multiply the cut-off frequency of the FIR filter
+        """
+        if x == 4:
+            return 1.004
+        elif x in (3, 1):
+            return 1.006
+        elif x in (-6, -4, -2, 2, 0):
+            return 1.01
+        else:
+            return 1.02
+
+    def decimation(self, x):
+        """
+        Required decimation for each filter
+        """
+        if x > 1:
+            return [1]
+        elif x > -2:
+            return [4]
+        elif x > -4:
+            return [4, 4]
+        elif x > -6:
+            return [4, 4, 4]
+        elif x == -6:
+            return [4, 4, 4, 4]
+        assert False, 'Overlooked decimation'
+
+
+class ThirdOctaveBankDesigner(FilterBankDesigner):
+
+    def __init__(self):
+        self.xs = list(range(-16, 14))
+        # Text corresponding to the nominal frequency
+        self.nominal_txt = ['25', '31.5', '40',
+                            '50', '63', '80',
+                            '100', '125', '160',
+                            '200', '250', '315',
+                            '400', '500', '630',
+                            '800', '1k', '1.25k',
+                            '1.6k', '2k', '2.5k',
+                            '3.15k', '4k', '5k',
+                            '6.3k', '8k', '10k',
+                            '12.5k', '16k', '20k']
+
+    @property
+    def b(self):
+        # Band division factor, 3 for one-third octave bands
+        return 3
+
+    def nominal(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]
+
+    @staticmethod
+    def decimation(x):
+        if x > 5:
+            return [1]
+        elif x > -1:
+            return [4]
+        elif x > -7:
+            return [4, 4]
+        elif x > -13:
+            return [4, 4, 4]
+        elif x > -17:
+            return [4, 4, 4, 4]
+        assert False, 'Bug: overlooked decimation'
+
+    @staticmethod
+    def fac_l(x):
+        if x in (-13, -7, -1, 5, 11, 12, 13):
+            return .995
+        elif x in (-12, -6, 0, 6):
+            return .98
+        else:
+            return .99
+
+    @staticmethod
+    def fac_u(x):
+        if x in (-14, -13, -8, -7, -1, -2, 3, 4, 5, 10, 11, 12):
+            return 1.005
+        elif x in (12, 13):
+            return 1.003
+        elif x in (12, -6, 0):
+            return 1.015
+        else:
+            return 1.01

lasp/filter/bandpass_limits.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE
+
+Description: Limit lines for class 1 octave band filter limits according to
+the ICS 17.140.50 standard.
+"""
+__all__ = ['G', 'fr', 'third_octave_band_limits', 'octave_band_limits']
+import numpy as np
+
+# Reference frequency
+fr = 1000.
+G = 10**(3/10)
+
+
+def third_octave_band_limits(x):
+    """
+    Returns the class 1 third octave band filter limits for filter designator
+    x.
+
+    Args:
+        x: Filter offset power from the reference frequency of 1000 Hz.
+
+    Returns:
+        freq, ulim, llim: Tuple of Numpy arrays containing the frequencyies,
+        upper and lower limits of the filter.
+    """
+    b = 3
+
+    fm = G**(x/b)*fr
+    plusinf = 20
+    f_ratio_pos = [1., 1.02667, 1.05575, 1.08746, 1.12202, 1.12202,
+                   1.29437, 1.88173, 3.05365, 5.39195, plusinf]
+
+    f_ratio_neg = [0.97402, 0.94719, 0.91958, 0.89125, 0.89125,
+                   0.77257, 0.53143, 0.32748, 0.18546, 1/plusinf]
+    f_ratio_neg.reverse()
+
+    f_ratio = f_ratio_neg + f_ratio_pos
+
+    mininf = -1e300
+
+    upper_limits_pos = [.3]*5 + [-2, -17.5, -42, -61, -70, -70]
+    upper_limits_neg = upper_limits_pos[:]
+    upper_limits_neg.reverse()
+    upper_limits = np.array(upper_limits_neg[:-1] + upper_limits_pos)
+
+    lower_limits_pos = [-.3, -.4, -.6, -1.3, -5, -5, mininf, mininf,
+                        mininf, mininf, mininf]
+    lower_limits_neg = lower_limits_pos[:]
+    lower_limits_neg.reverse()
+    lower_limits = np.array(lower_limits_neg[:-1] + lower_limits_pos)
+
+    freqs = fm*np.array(f_ratio)
+
+    return freqs, upper_limits, lower_limits
+
+
+def octave_band_limits(x):
+
+    b = 1
+
+    # Exact midband frequency
+    fm = G**(x/b)*fr
+
+    G_power_values_pos = [0, 1/8, 1/4, 3/8, 1/2, 1/2, 1, 2, 3, 4]
+    G_power_values_neg = [-i for i in G_power_values_pos]
+    G_power_values_neg.reverse()
+    G_power_values = G_power_values_neg[:-1] + G_power_values_pos
+
+    mininf = -1e300
+
+    lower_limits_pos = [-0.3, -0.4, -0.6, -1.3, -5.0, -5.0] + 4*[mininf]
+    lower_limits_neg = lower_limits_pos[:]
+    lower_limits_neg.reverse()
+    lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos)
+
+    upper_limits_pos = [0.3]*5 + [-2, -17.5, -42, -61, -70]
+    upper_limits_neg = upper_limits_pos[:]
+    upper_limits_neg.reverse()
+    upper_limits = np.asarray(upper_limits_neg[:-1] + upper_limits_pos)
+
+    freqs = fm*G**np.asarray(G_power_values)
+
+    return freqs, upper_limits, lower_limits
+
+
+if __name__ == '__main__':
+
+    from asceefigs.plot import close, Figure
+    close('all')
+    freqs, upper_limits, lower_limits = octave_band_limits(0)
+
+    f = Figure()
+    f.semilogx(freqs, lower_limits)
+    f.semilogx(freqs, upper_limits)
+
+    f.ylim(-80, 1)

lasp/filter/dec_fir.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE
+
+Description: Decimation filter design.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from fir_design import freqResponse, lowpass_fir_design
+
+L = 128  # Filter order
+fs = 48000.  # Sampling frequency
+
+d = 4 # Decimation factor
+
+fd = fs/d  # Decimated sampling frequency
+fc = fd/2/1.5  # Filter cut-off frequency
+
+fir = lowpass_fir_design(L, fs, fc)
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+
+freq = np.logspace(1, np.log10(fs/2), 1000)
+
+H = freqResponse(fs, freq, fir)
+dBH = 20*np.log10(np.abs(H))
+ax.plot(freq, dBH)
+ax.axvline(fs/2, color='green')
+ax.axvline(fd/2, color='red')
+
+# Index of Nyquist freq
+fn_index = np.where(freq <= fd/2)[0][-1]
+
+dBHmax_above_Nyq = np.max(dBH[fn_index:])
+
+print(f"Reduction above Nyquist: {dBHmax_above_Nyq} dB")
+plt.show()

lasp/filter/fir_design.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE
+
+Description: Designs octave band FIR filters from 16Hz to 16 kHz for a sampling
+frequency of 48 kHz.
+"""
+# from asceefigs.plot import Bode, close, Figure
+__all__ = ['freqResponse', 'bandpass_fir_design', 'lowpass_fir_design',
+           'arbitrary_fir_design']
+
+import numpy as np
+from scipy.signal import freqz, hann, firwin2
+
+
+def freqResponse(fs, freq, coefs_b, coefs_a=1.):
+    """
+    Computes the frequency response of the filter defined with the filter
+    coefficients:
+
+    Args:
+        fs: Sampling frequency [Hz]
+        freq: Array of frequencies to compute the response for
+        coefs_b: Forward coefficients (FIR coefficients)
+        coefs_a: Feedback coefficients (IIR)
+
+    Returns:
+        Complex frequency response for frequencies given in array
+    """
+    Omg = 2*np.pi*freq/fs
+
+    w, H = freqz(coefs_b, coefs_a, worN=Omg)
+    return H
+
+
+def bandpass_fir_design(L, fs, fl, fu, window=hann):
+    """
+    Construct a bandpass filter
+    """
+    assert fs/2 > fu, "Nyquist frequency needs to be higher than upper cut-off"
+    assert fu > fl, "Cut-off needs to be lower than Nyquist freq"
+
+    Omg2 = 2*np.pi*fu/fs
+    Omg1 = 2*np.pi*fl/fs
+
+    fir = np.empty(L, dtype=float)
+
+    # First Create ideal band-pass filter
+    fir[L//2] = (Omg2-Omg1)/np.pi
+
+    for n in range(1, L//2):
+        fir[n+L//2] = (np.sin(n*Omg2)-np.sin(n*Omg1))/(n*np.pi)
+        fir[L//2-n] = (np.sin(n*Omg2)-np.sin(n*Omg1))/(n*np.pi)
+
+    win = window(L, True)
+    fir_win = fir*win
+
+    return fir_win
+
+
+def lowpass_fir_design(L, fs, fc, window=hann):
+    assert fs/2 > fc, "Nyquist frequency needs to be higher" \
+                      " than upper cut-off"
+
+    Omgc = 2*np.pi*fc/fs
+    fir = np.empty(L, dtype=float)
+
+    # First Create ideal band-pass filter
+    fir[L//2] = Omgc/np.pi
+
+    for n in range(1, L//2):
+        fir[n+L//2] = np.sin(n*Omgc)/(n*np.pi)
+        fir[L//2-n] = np.sin(n*Omgc)/(n*np.pi)
+
+    win = window(L, True)
+    fir_win = fir*win
+
+    return fir_win
+
+
+def arbitrary_fir_design(fs, L, freq, amps, window='hann'):
+    """
+    Last frequency of freq should be fs/2
+    """
+    return firwin2(L, freq, amps, fs=fs, window=window)

lasp/filter/freqweighting_fir.py

@@ -6,11 +6,10 @@ Author: J.A. de Jong - ASCEE
 Description: Filter design for frequency weightings A and C.
 """
 from .fir_design import freqResponse, arbitrary_fir_design
-from ..lasp_common import FreqWeighting
 import numpy as np
 
-__all__ = ['A','C','A_fir_design','C_fir_design',
+__all__ = ['A', 'C', 'A_fir_design', 'C_fir_design',
-           'show_Afir','show_Cfir']
+           'show_Afir', 'show_Cfir']
 
 fr = 1000.
 fL = 10**1.5
@@ -53,7 +52,7 @@ def C_uncor(f):
     Computes the uncorrected frequency response of the C-filter
     """
     fsq = f**2
-    num  = f4sq*fsq
+    num = f4sq*fsq
     denom1 = (fsq+f1**2)
     denom2 = (fsq+f4**2)
     return num/(denom1*denom2)
@@ -69,79 +68,79 @@ def C(f):
 
 def A_fir_design():
     fs = 48000.
-    freq_design = np.linspace(0,17e3,3000)
+    freq_design = np.linspace(0, 17e3, 3000)
     freq_design[-1] = fs/2
-    amp_design  = A(freq_design)
+    amp_design = A(freq_design)
     amp_design[-1] = 0.
 
-    L = 2048 # Filter order
+    L = 2048  # Filter order
-    fir = arbitrary_fir_design(fs,L,freq_design,amp_design,
+    fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
                                window='rectangular')
     return fir
 
 
 def C_fir_design():
     fs = 48000.
-    freq_design = np.linspace(0,17e3,3000)
+    freq_design = np.linspace(0, 17e3, 3000)
     freq_design[-1] = fs/2
-    amp_design  = C(freq_design)
+    amp_design = C(freq_design)
     amp_design[-1] = 0.
 
-    L = 2048 # Filter order
+    L = 2048  # Filter order
-    fir = arbitrary_fir_design(fs,L,freq_design,amp_design,
+    fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
                                window='rectangular')
     return fir
 
 
 def show_Afir():
-    from asceefigs.plot import Bode, close, Figure
+    from asceefigs.plot import close, Figure
     close('all')
 
     fs = 48000.
-    freq_design = np.linspace(0,17e3,3000)
+    freq_design = np.linspace(0, 17e3, 3000)
     freq_design[-1] = fs/2
-    amp_design  = A(freq_design)
+    amp_design = A(freq_design)
     amp_design[-1] = 0.
     firs = []
 
-#    firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hamming'))
+    # firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hamming'))
-#    firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hann'))
+    # firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hann'))
     firs.append(A_fir_design())
-    #from scipy.signal import iirdesign
+    # from scipy.signal import iirdesign
-    #b,a = iirdesign()
+    # b,a = iirdesign()
-    freq_check = np.logspace(0,np.log10(fs/2),5000)
+    freq_check = np.logspace(0, np.log10(fs/2), 5000)
-    f=Figure()
+    f = Figure()
 
-    f.semilogx(freq_check,20*np.log10(A(freq_check)))
+    f.semilogx(freq_check, 20*np.log10(A(freq_check)))
     for fir in firs:
-        H = freqResponse(fs,freq_check, fir)
+        H = freqResponse(fs, freq_check, fir)
-        f.plot(freq_check,20*np.log10(np.abs(H)))
+        f.plot(freq_check, 20*np.log10(np.abs(H)))
 
-    f.fig.get_axes()[0].set_ylim(-75,3)
+    f.fig.get_axes()[0].set_ylim(-75, 3)
 
 
 def show_Cfir():
-    from asceefigs.plot import Bode, close, Figure
+    from asceefigs.plot import close, Figure
     close('all')
 
     fs = 48000.
-    freq_design = np.linspace(0,17e3,3000)
+    freq_design = np.linspace(0, 17e3, 3000)
     freq_design[-1] = fs/2
-    amp_design  = C(freq_design)
+    amp_design = C(freq_design)
     amp_design[-1] = 0.
     firs = []
 
-#    firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hamming'))
+    # firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hamming'))
-#    firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hann'))
+    # firs.append(arbitrary_fir_design(fs,L,freq_design,amp_design,window='hann'))
     firs.append(C_fir_design())
-    #from scipy.signal import iirdesign
+    # from scipy.signal import iirdesign
-    #b,a = iirdesign()
+    # b,a = iirdesign()
-    freq_check = np.logspace(0,np.log10(fs/2),5000)
+    freq_check = np.logspace(0, np.log10(fs/2), 5000)
-    f=Figure()
+    f = Figure()
 
-    f.semilogx(freq_check,20*np.log10(C(freq_check)))
+    f.semilogx(freq_check, 20*np.log10(C(freq_check)))
     for fir in firs:
-        H = freqResponse(fs,freq_check, fir)
+        H = freqResponse(fs, freq_check, fir)
-        f.plot(freq_check,20*np.log10(np.abs(H)))
+        f.plot(freq_check, 20*np.log10(np.abs(H)))
 
-    f.fig.get_axes()[0].set_ylim(-30,1)
+    f.fig.get_axes()[0].set_ylim(-30, 1)

lasp/filter_design/fir_design.py

@@ -1,70 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""!
-Author: J.A. de Jong - ASCEE
-
-Description: Designs octave band FIR filters from 16Hz to 16 kHz for a sampling
-frequency of 48 kHz.
-"""
-#from asceefigs.plot import Bode, close, Figure
-
-import numpy as np
-from scipy.signal import freqz, hann, firwin2
-from matplotlib.pyplot import figure, close
-
-def freqResponse(fs, freq, fir_coefs_b, fir_coefs_a=1.):
-    """!
-    Computes the frequency response of the filter defined with filter_coefs
-    """
-    Omg = 2*np.pi*freq/fs
-
-    w, H = freqz(fir_coefs_b,fir_coefs_a,worN = Omg)
-    return H
-
-def bandpass_fir_design(L,fs,fl,fu, window = hann):
-    """
-    Construct a bandpass filter
-    """
-    assert fs/2 > fu, "Nyquist frequency needs to be higher than upper cut-off"
-    assert fu > fl, "Cut-off needs to be lower than Nyquist freq"
-    
-    Omg2 = 2*np.pi*fu/fs
-    Omg1 = 2*np.pi*fl/fs
-        
-    fir = np.empty(L, dtype=float)
-
-    # First Create ideal band-pass filter
-    fir[L//2] = (Omg2-Omg1)/np.pi
-        
-    for n in range(1,L//2):
-        fir[n+L//2] = (np.sin(n*Omg2)-np.sin(n*Omg1))/(n*np.pi)
-        fir[L//2-n] = (np.sin(n*Omg2)-np.sin(n*Omg1))/(n*np.pi)
-        
-    win = window(L,True)
-    fir_win = fir*win
-        
-    return fir_win    
-
-def lowpass_fir_design(L,fs,fc,window = hann):
-    assert fs/2 > fc, "Nyquist frequency needs to be higher than upper cut-off"
-    
-    Omgc = 2*np.pi*fc/fs
-    fir = np.empty(L, dtype=float)
-
-    # First Create ideal band-pass filter
-    fir[L//2] = Omgc/np.pi
-        
-    for n in range(1,L//2):
-        fir[n+L//2] = np.sin(n*Omgc)/(n*np.pi)
-        fir[L//2-n] = np.sin(n*Omgc)/(n*np.pi)
-        
-    win = window(L,True)
-    fir_win = fir*win
-        
-    return fir_win    
-
-def arbitrary_fir_design(fs,L,freq,amps,window='hann'):
-    """
-    Last frequency of freq should be fs/2
-    """
-    return firwin2(L,freq,amps,fs=fs,window=window)

lasp/filter_design/octave_filter_design_fir.py

@@ -1,183 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""!
-Author: J.A. de Jong - ASCEE
-
-Description: Designs FIR octave band FIR filters from 16Hz to 16 kHz for a
-sampling frequency of 48 kHz.
-"""
-import numpy as np
-from octave_filter_limits import octave_band_limits, G, fr
-from matplotlib.pyplot import figure, close
-from fir_design import freqResponse, bandpass_fir_design
-
-L = 256 # Filter order
-fs = 48000. # Sampling frequency
-
-showfig = False
-#showfig = True
-maxdec = 3
-close('all')
-
-filters = {}
-decimation = {}
-b = 1
-# Text corresponding to the nominal frequency
-nominals = {4:'16k',
-            3:'8k',
-            2:'4k',
-            1:'2k',
-            0:'1k',
-            -1:'500',
-            -2:'250',
-            -3:'125',
-            -4:'63',
-            -5:'31.5',
-            -6:'16'}
-
-# Factor with which to multiply the cut-on frequency of the FIR filter
-cut_fac_l = {
-4:0.995,
-3:0.99,
-2:0.98,
-1:0.99,
-0:0.98,
--1:0.96,
--2:.99,
--3:0.96,
--4:0.96,
--5:0.96,
--6:0.96
-}
-
-# Factor with which to multiply the cut-off frequency of the FIR filter
-cut_fac_u = {
-4:1.004,
-3:1.006,
-2:1.01,
-1:1.006,
-0:1.01,
--1:1.02,
--2:1.006,
--3:1.02,
--4:1.02,
--5:1.02,
--6:1.02
-}
-
-# Required decimation for each filter
-decimation = {
-4:[1],
-3:[1],
-2:[1],
-1:[4],
-0:[4],
--1:[4],
--2:[4,4],
--3:[4,4],
--4:[4,4,4],
--5:[4,4,4],
--6:[4,4,4,4]
-        }
-
-# Generate the header file
-with open('../c/lasp_octave_fir.h','w') as hfile:
-    hfile.write("""// lasp_octave_fir.h
-//
-// Author: J.A. de Jong - ASCEE
-//
-// Description: This is file is automatically generated from the
-// octave_filter_design.py Python file.
-//////////////////////////////////////////////////////////////////////
-#pragma once
-#ifndef LASP_OCTAVE_FIR_H
-#define LASP_OCTAVE_FIR_H
-#include "lasp_types.h"
-
-#define MAXDEC (%i) /// Size of the decimation factor array
-#define OCTAVE_FIR_LEN (%i) /// Filter length
-
-typedef struct {
-    const char* nominal; /// Pointer to the nominal frequency text
-    int x; /// 1000*G^x is the exact midband frequency, where G = 10^(3/10)
-    int decimation_fac[MAXDEC]; // Array with decimation factors that need to be
-                                // applied prior to filtering
-    d h[OCTAVE_FIR_LEN];
-} OctaveFIR;
-
-extern __thread OctaveFIR OctaveFIRs[%i];
-        
-#endif // LASP_OCTAVE_FIR_H
-//////////////////////////////////////////////////////////////////////
-""" %(maxdec,L,len(decimation)))
-    
-# Generate the source code file
-with open('../c/lasp_octave_fir.c','w') as cfile:
-    cfile.write("""// octave_fir.c
-//
-// Author: J.A. de Jong - ASCEE
-// 
-// Description:
-// This is an automatically generated file containing filter coefficients
-// for octave filter banks.
-//////////////////////////////////////////////////////////////////////
-#include "lasp_octave_fir.h"
-
-__thread OctaveFIR OctaveFIRs[%i] = {
-""" %(len(decimation)))
-    struct = ''
-    for x in range(4,-7,-1):
-        
-        if showfig:
-            fig = figure()
-            ax = fig.add_subplot(111)
-        
-        dec = decimation[x]
-        d = np.prod(dec)
-        
-        struct += '\n    {"%s",%i, {' %(nominals[x],x)
-        for i in range(maxdec):
-            try:
-                struct += "%i," % dec[i]
-            except:
-                struct += "0,"
-        struct = struct[:-1] # Strip off last,
-        struct += '},'
-        
-        fd = fs/d
-        
-        # Exact midband frequency
-        fm = G**(x)*fr
-        
-        # Cut-on frequency of the filter
-        f1 = fm*G**(-1/(2*b))*cut_fac_l[x]
-        # 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 += "{ "
-        for i in fir:
-            struct += "%0.16e," %i
-        struct = struct[:-1] + " }},"
-        
-        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)))
-            
-            freq,ulim,llim = octave_band_limits(x)
-            
-            ax.semilogx(freq,ulim)
-            ax.semilogx(freq,llim)
-            ax.set_title('x = %i, fnom = %s' %(x,nominals[x]) )
-        
-            ax.set_ylim(-10,1)
-            ax.set_xlim(f1/1.1,fd)
-            ax.axvline(fd/2)
-            ax.axvline(fc,color='red')
-    struct+="\n};"
-    cfile.write(struct)    
-    

lasp/filter_design/octave_filter_limits.py

@@ -1,52 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""!
-Author: J.A. de Jong - ASCEE
-
-Description: Limit lines for class 1 octave band filter limits according to
-the ICS 17.140.50 standard.
-"""
-__all__ = ['G','fr','octave_band_limits']
-import numpy as np
-
-# Reference frequency
-fr = 1000.
-G = 10**(3/10)
-
-def octave_band_limits(x):
-    
-    # Exact midband frequency
-    fm = G**(x)*fr
-
-    G_power_values_pos = [0,1/8,1/4,3/8,1/2,1/2,1,2,3,4]
-    G_power_values_neg = [-i for i in G_power_values_pos]
-    G_power_values_neg.reverse()
-    G_power_values = G_power_values_neg[:-1] + G_power_values_pos
-    
-    mininf = -1e300
-    
-    lower_limits_pos = [-0.3,-0.4,-0.6,-1.3,-5.0,-5.0]+ 4*[mininf]
-    lower_limits_neg =  lower_limits_pos[:]
-    lower_limits_neg.reverse()
-    lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos)
-    
-    upper_limits_pos = [0.3]*5 + [-2,-17.5,-42,-61,-70]
-    upper_limits_neg =  upper_limits_pos[:]
-    upper_limits_neg.reverse()
-    upper_limits = np.asarray(upper_limits_neg[:-1] + upper_limits_pos)
-    
-    freqs = fm*G**np.asarray(G_power_values)
-
-    return freqs,upper_limits,lower_limits
-
-if __name__ == '__main__':
-    
-    from asceefigs.plot import close, Figure
-    close('all')
-    freqs,upper_limits,lower_limits = octave_band_limits(0)
-
-    f = Figure()
-    f.semilogx(freqs,lower_limits)
-    f.semilogx(freqs,upper_limits)
-    
-    f.ylim(-80,1)

lasp/lasp_octavefilter.py

@@ -0,0 +1,154 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE
+
+Provides the FIR implementation of the octave filter bank
+
+"""
+__all__ = ['OctaveFilterBank']
+from .filter.bandpass_fir import OctaveBankDesigner, ThirdOctaveBankDesigner
+from .wrappers import Decimator, FilterBank as pyxFilterBank
+import numpy as np
+
+
+class FilterBank:
+    """
+    Single channel octave filter bank implementation
+    """
+
+    def __init__(self, fs, designer):
+        """
+        Initialize a OctaveFilterBank object.
+
+        Args:
+            fs: Sampling frequency of base signal
+            designer: FIR Filter designer for filterbank
+
+        """
+        assert np.isclose(fs, 48000), "Only sampling frequency" \
+            " available is 48 kHz"
+
+        self.fs = fs
+
+        maxdecimation = designer.decimation(designer.xs[0])
+        self.decimators = []
+        for dec in maxdecimation:
+            self.decimators.append(Decimator(1, dec))
+
+        xs_d1 = []
+        xs_d4 = []
+        xs_d16 = []
+        xs_d64 = []
+        xs_d256 = []
+
+        self.filterbanks = []
+        # Sort the x values in categories according to the required decimation
+        for x in designer.xs:
+            dec = designer.decimation(x)
+            if len(dec) == 1 and dec[0] == 1:
+                xs_d1.append(x)
+            elif len(dec) == 1 and dec[0] == 4:
+                xs_d4.append(x)
+            elif len(dec) == 2:
+                xs_d16.append(x)
+            elif len(dec) == 3:
+                xs_d64.append(x)
+            elif len(dec) == 4:
+                xs_d256.append(x)
+            else:
+                raise ValueError(f'No decimation found for x={x}')
+
+        xs_all = [xs_d1, xs_d4, xs_d16, xs_d64, xs_d256]
+        for xs in xs_all:
+            nominals = []
+            firs = np.empty((designer.L, len(xs)), order='F')
+            for i, x in enumerate(xs):
+                #  These are the filters that do not require lasp_decimation
+                #  prior to filtering
+                nominals.append(designer.nominal(x))
+                firs[:, i] = designer.createFilter(fs, x)
+            filterbank = {'fb': pyxFilterBank(firs, 1024),
+                          'xs': xs,
+                          'nominals': nominals}
+            self.filterbanks.append(filterbank)
+
+        # Sample input counter.
+        self.N = 0
+
+        # Filter output counters
+        self.dec = [1, 4, 16, 64, 256]
+        Pdec = 128
+        Pfil = 256
+
+        # Initial filtered number
+        self.Nf = [-Pfil, -(Pdec/4 + Pfil),
+                   -(Pfil + Pdec/4 + Pdec/16),
+                   -(Pfil + Pdec/64 + Pdec/4 + Pdec/16),
+                   -(Pfil + Pdec/256 + Pdec/64 + Pdec/4 + Pdec/16)]
+
+        self.delay_set = [False, False, False, False, False]
+
+    def filterd(self, dec_stage, data):
+        """
+        Filter data for a given decimation stage
+
+        Args:
+            dec_stage: decimation stage
+            data: Pre-filtered data
+        """
+        output = {}
+        if data.shape[0] == 0:
+            return output
+
+        filtered = self.filterbanks[dec_stage]['fb'].filter_(data)
+        Nf = filtered.shape[0]
+        if Nf > 0:
+            dec = self.dec[dec_stage]
+            fd = self.fs/dec
+
+            oldNf = self.Nf[dec_stage]
+            tstart = oldNf/fd
+            tend = tstart + Nf/fd
+            t = np.linspace(tstart, tend, Nf, endpoint=False)
+            self.Nf[dec_stage] += Nf
+            for i, nom in enumerate(self.filterbanks[dec_stage]['nominals']):
+                output[nom] = {'t': t, 'data': filtered[:, i]}
+        return output
+
+    def filter_(self, data):
+        """
+        Filter input data
+        """
+        assert data.ndim == 2
+        assert data.shape[1] == 1, "invalid number of channels, should be 1"
+
+        if data.shape[0] == 0:
+            return {}
+
+        # Output given as a dictionary with x as the key
+        output = {}
+        self.N += data.shape[0]
+
+        output = {**output, **self.filterd(0, data)}
+
+        for i in range(len(self.decimators)):
+            dec_stage = i+1
+            if data.shape[0] > 0:
+                # Apply a decimation stage
+                data = self.decimators[i].decimate(data)
+                output = {**output, **self.filterd(dec_stage, data)}
+
+        return output
+
+
+class OctaveFilterBank(FilterBank):
+    def __init__(self, fs):
+        designer = OctaveBankDesigner()
+        super().__init__(fs, designer)
+
+
+class ThirdOctaveFilterBank(FilterBank):
+    def __init__(self, fs):
+        designer = ThirdOctaveBankDesigner()
+        super().__init__(fs, designer)

lasp/wrappers.pyx

@@ -304,9 +304,10 @@ cdef class FilterBank:
         if self.fb:
             FilterBank_free(self.fb)
 
-    def filter_(self,d[:] input_):
+    def filter_(self,d[::1, :] input_):
+        assert input_.shape[1] == 1
         cdef dmat input_vd = dmat_foreign_data(input_.shape[0],1,
-                                             &input_[0],False)
+                                             &input_[0, 0],False)
 
         cdef dmat output = FilterBank_filter(self.fb,&input_vd)
 
@@ -340,6 +341,9 @@ cdef class Decimator:
 
     def decimate(self,d[::1,:] samples):
         assert samples.shape[1] == self.nchannels,'Invalid number of channels'
+        if samples.shape[0] == 0:
+            return np.zeros((0, self.nchannels))
+
         cdef dmat d_samples = dmat_foreign_data(samples.shape[0],
                                                 samples.shape[1],
                                                 &samples[0,0],

test/bandpass_test_1.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE
+
+"""
+import numpy as np
+from lasp.filter.bandpass_limits import (third_octave_band_limits,
+                                         octave_band_limits, G, fr)
+
+from lasp.filter.bandpass_fir import ThirdOctaveBankDesigner, \
+                                     OctaveBankDesigner
+import matplotlib.pyplot as plt
+
+# Adjust these settings
+b = 1  # or three
+zoom = False  # or True
+
+if b == 3:
+    bands = ThirdOctaveBankDesigner()
+elif b == 1:
+    bands = OctaveBankDesigner()
+else:
+    raise ValueError('b should be 1 or 3')
+
+
+for x in bands.xs:
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+    fs = 48000.
+    dec = np.prod(bands.decimation(x))
+    fd = fs/dec
+    fc = fd/2/1.4
+
+    freq = np.logspace(np.log10(1), np.log10(fd), 5000)
+    H = bands.freqResponse(fs, x, freq)
+    dBH = 20*np.log10(np.abs(H))
+    ax.semilogx(freq, dBH)
+
+    if b == 1:
+        freq, ulim, llim = octave_band_limits(x)
+    else:
+        freq, ulim, llim = third_octave_band_limits(x)
+
+    ax.semilogx(freq, llim)
+    ax.semilogx(freq, ulim)
+    ax.set_title(f'x = {x}, fnom = {bands.nominal(x)}')
+
+    if zoom:
+        ax.set_xlim(bands.fl(x)/1.1, bands.fu(x)*1.1)
+        ax.set_ylim(-15, 1)
+    else:
+        ax.set_ylim(-75, 1)
+        ax.set_xlim(10, fd)
+
+    ax.axvline(fd/2)
+    if dec > 1:
+        ax.axvline(fc, color='red')
+
+    ax.legend(['Filter frequency response',
+               'Lower limit from standard',
+               'Upper limit from standard',
+               'Nyquist frequency after decimation',
+               'Decimation filter cut-off frequency'], fontsize=8)
+
+plt.show()

test/test_bars.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import sys
+from lasp.plot import BarScene
+from PySide import QtGui
+from PySide.QtCore import QTimer
+from lasp.lasp_gui_tools import Branding, ASCEEColors
+import numpy as np
+import PySide.QtOpenGL as gl
+
+
+def main():
+    app = QtGui.QApplication(sys.argv)  # A new instance of QApplication
+    app.setFont(Branding.font())
+    pix = QtGui.QPixmap(':img/img/lasp_logo_640.png')
+    splash = QtGui.QSplashScreen(pixmap=pix)
+    splash.show()
+    mw = QtGui.QGraphicsView()
+    glwidget = gl.QGLWidget()
+    mw.setViewport(glwidget)
+    bs = BarScene(None, np.array([10, 20, 300]), 2, ylim=(0, 1))
+    mw.setScene(bs)
+
+    bs.set_ydata(np.array([[.1, .2],
+                           [.7, .8],
+                           [.9, 1]]))
+
+    # timer = QTimer.
+    print(ASCEEColors.bggreen.getRgb())
+    mw.show()                         # Show the form
+    splash.finish(mw)
+    app.exec_()                         # and execute the app
+
+
+if __name__ == '__main__':
+    main()                    # run the main function