Cleanup fir filter design code. Removed Beamforming reference i test_fft.py

.gitignore

@@ -16,4 +16,5 @@ test/test_fft
 test/test_math
 test/test_workers
 doc
-
+LASP.egg-info
+lasp_octave_fir.*

lasp/c/CMakeLists.txt

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

lasp/fir_design/fir_design.py

@@ -0,0 +1,65 @@
+#!/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
+from matplotlib.pyplot import figure, close
+
+def freqResponse(fir_coefs, freq, fs):
+    """!
+    Computes the frequency response of the filter defined with filter_coefs
+    """
+    Omg = 2*np.pi*freq/fs
+
+    w, H = freqz(fir_coefs,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    
+    

lasp/fir_design/octave_filter_design_fir.py

@@ -0,0 +1,180 @@
+#!/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.98,
+-4:0.96,
+-5:0.98,
+-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.01,
+-4:1.02,
+-5:1.01,
+-6:1.02
+}
+
+# Required decimation for each filter
+decimation = {
+4:[1],
+3:[1],
+2:[1],
+1:[4],
+0:[4],
+-1:[4],
+-2:[4,8],
+-3:[4,8],
+-4:[4,8],
+-5:[4,8,4],
+-6:[4,8,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]
+
+        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(f2*3),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,f2*1.1)
+    
+    struct+="\n};"
+    cfile.write(struct)    
+    

lasp/fir_design/octave_filter_limits.py

@@ -0,0 +1,52 @@
+#!/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)

test/fft_test.py

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