Working on SOS implementation of frequency weighting filters

lasp/filter/soundpressureweighting.py

@@ -7,103 +7,137 @@ Description: Filter design for frequency weighting curves (i.e. A and C
 weighting)
 """
 from .fir_design import freqResponse, arbitrary_fir_design
+from scipy.signal import bilinear_zpk, zpk2sos
 import numpy as np
 
-__all__ = ['A', 'C', 'A_fir_design', 'C_fir_design',
+__all__ = ['SPLFilterDesigner']
-           'show_Afir', 'show_Cfir']
-
-fr = 1000.
-fL = 10**1.5
-fH = 10**3.9
-fLsq = fL**2
-fHsq = fH**2
-frsq = fr**2
-fA = 10**2.45
-D = np.sqrt(.5)
-
-b = (1/(1-D))*(frsq+fLsq*fHsq/frsq-D*(fLsq+fHsq))
-c = fLsq*fHsq
-f2 = (3-np.sqrt(5.))/2*fA
-f3 = (3+np.sqrt(5.))/2*fA
-f1 = np.sqrt((-b-np.sqrt(b**2-4*c))/2)
-f4 = np.sqrt((-b+np.sqrt(b**2-4*c))/2)
-f4sq = f4**2
 
 
-def A_uncor(f):
+class SPLFilterDesigner:
-    """
+    fr = 1000.
-    Computes the uncorrected frequency response of the A-filter
+    fL = 10**1.5
-    """
+    fH = 10**3.9
-    fsq = f**2
-    num = f4sq*fsq**2
-    denom1 = (fsq+f1**2)
-    denom2 = np.sqrt((fsq+f2**2)*(fsq+f3**2))*(fsq+f4sq)
-
-    return (num/(denom1*denom2))
 
 
-def A(f):
+    fLsq = fL**2
-    """
+    fHsq = fH**2
-    Computes the linear A-weighting freqency response. Hence, to obtain
+    frsq = fr**2
-    A-weighted values, the *amplitude* need to be multiplied with this value.
+    fA = 10**2.45
-    Hence, to correct dB levels, the value of 20*log(A) needs to be added to
+    D = np.sqrt(.5)
-    the level
 
-    Args:
+    b = (1/(1-D))*(frsq+fLsq*fHsq/frsq-D*(fLsq+fHsq))
-        f: Frequency array to compute values for
+    c = fLsq*fHsq
-    Returns:
+    f2 = (3-np.sqrt(5.))/2*fA
-        A(f) for each frequency
+    f3 = (3+np.sqrt(5.))/2*fA
-    """
+    f1 = np.sqrt((-b-np.sqrt(b**2-4*c))/2)
-    Auncor = A_uncor(f)
+    f4 = np.sqrt((-b+np.sqrt(b**2-4*c))/2)
-    A1000 = A_uncor(1000.)
+    f4sq = f4**2
-    return Auncor/A1000
+
+    def _A_uncor(self, f):
+        """
+        Computes the uncorrected frequency response of the A-filter
+
+        Args:
+            f: Frequency (array, float)
+
+        Returns:
+            Linear filter transfer function
+        """
+        fsq = f**2
+        num = self.f4sq*fsq**2
+        denom1 = (fsq+self.f1**2)
+        denom2 = np.sqrt((fsq+self.f2**2)*(fsq+self.f3**2))*(fsq+self.f4sq)
+
+        return (num/(denom1*denom2))
 
 
-def C_uncor(f):
+    def A(self, f):
-    """
+        """
-    Computes the uncorrected frequency response of the C-filter
+        Computes the linear A-weighting freqency response. Hence, to obtain
-    """
+        A-weighted values, the *amplitude* need to be multiplied with this value.
-    fsq = f**2
+        Hence, to correct dB levels, the value of 20*log(A) needs to be added to
-    num = f4sq*fsq
+        the level
-    denom1 = (fsq+f1**2)
+
-    denom2 = (fsq+f4**2)
+        Args:
-    return num/(denom1*denom2)
+            f: Frequency array to compute values for
+        Returns:
+            A(f) for each frequency
+        """
+        Auncor = self._A_uncor(f)
+        A1000 = self._A_uncor(self.fr)
+        return Auncor/A1000
 
 
-def C(f):
+    def _C_uncor(self, f):
-    """
+        """
-    Computes the linear A-weighting freqency response
+        Computes the uncorrected frequency response of the C-filter
-    """
+        """
-    Cuncor = C_uncor(f)
+        fsq = f**2
-    C1000 = C_uncor(1000.)
+        num = self.f4sq*fsq
-    return Cuncor/C1000
+        denom1 = (fsq+self.f1**2)
+        denom2 = (fsq+self.f4**2)
+        return num/(denom1*denom2)
 
 
-def A_fir_design():
+    def C(self, f):
-    fs = 48000.
+        """
-    freq_design = np.linspace(0, 17e3, 3000)
+        Computes the linear A-weighting freqency response
-    freq_design[-1] = fs/2
+        """
-    amp_design = A(freq_design)
+        Cuncor = self._C_uncor(f)
-    amp_design[-1] = 0.
+        C1000 = self._C_uncor(self.fr)
-
+        return Cuncor/C1000
-    L = 2048  # Filter order
-    fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
-                               window='rectangular')
-    return fir
 
 
-def C_fir_design():
+    def A_fir_design(self, fs):
-    fs = 48000.
-    freq_design = np.linspace(0, 17e3, 3000)
-    freq_design[-1] = fs/2
-    amp_design = C(freq_design)
-    amp_design[-1] = 0.
 
-    L = 2048  # Filter order
+        assert int(fs) == 48000
-    fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
+        freq_design = np.linspace(0, 17e3, 3000)
-                               window='rectangular')
+        freq_design[-1] = fs/2
-    return fir
+        amp_design = self.A(freq_design)
+        amp_design[-1] = 0.
 
+        L = 2048  # Filter order
+        fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
+                                   window='rectangular')
+        return fir
+
+
+    def C_fir_design(self, fs):
+        assert int(fs) == 48000
+        fs = 48000.
+        freq_design = np.linspace(0, 17e3, 3000)
+        freq_design[-1] = fs/2
+        amp_design = C(freq_design)
+        amp_design[-1] = 0.
+
+        L = 2048  # Filter order
+        fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
+                                   window='rectangular')
+        return fir
+
+    def C_Sos_design(self, fs):
+        """
+        Create filter coefficients of the C-weighting filter. Uses the bilinear
+        transform to convert the analog filter to a digital one.
+
+        Args:
+            fs: Sampling frequency [Hz]
+
+        Returns:
+            Sos: Second order sections
+        """
+        
+        p1 = 2*np.pi*self.f1
+        p4 = 2*np.pi*self.f4
+        zeros_analog = [0,0]
+        poles_analog = [p1, p1, p4, p4]
+        k_analog = p4**2/self._C_uncor(self.fr)
+
+        z, p, k = bilinear_zpk(zeros_analog, poles_analog, k_analog, fs)
+        sos = zpk2sos(z, p, k)
+        return sos
+        # return z, p, k
+        # return zeros_analog, poles_analog, k_analog
 
 def show_Afir():
     from asceefig.plot import Figure
@@ -155,3 +189,4 @@ def show_Cfir():
         f.plot(freq_check, 20*np.log10(np.abs(H)))
 
     f.fig.get_axes()[0].set_ylim(-30, 1)
+

lasp/lasp_weighcal.py

@@ -5,7 +5,7 @@ Weighting and calibration filter in one
 @author: J.A. de Jong - ASCEE
 """
 from .lasp_common import FreqWeighting
-from .filter import (A, C, arbitrary_fir_design, freqResponse as frp)
+from .filter import SPLFilterDesigner
 from lasp.lasp_config import ones, empty
 from .wrappers import FilterBank
 import numpy as np