2018-04-23 07:29:21 +00:00
|
|
|
#!/usr/bin/env python3
|
|
|
|
# -*- coding: utf-8 -*-
|
|
|
|
"""!
|
|
|
|
Author: J.A. de Jong - ASCEE
|
|
|
|
|
2020-01-03 21:00:50 +00:00
|
|
|
Description: Filter design for frequency weighting curves (i.e. A and C
|
|
|
|
weighting)
|
2018-04-23 07:29:21 +00:00
|
|
|
"""
|
|
|
|
from .fir_design import freqResponse, arbitrary_fir_design
|
2020-01-17 21:10:42 +00:00
|
|
|
from scipy.signal import bilinear_zpk, zpk2sos
|
2018-04-23 07:29:21 +00:00
|
|
|
import numpy as np
|
|
|
|
|
2020-01-17 21:10:42 +00:00
|
|
|
__all__ = ['SPLFilterDesigner']
|
|
|
|
|
|
|
|
|
|
|
|
class SPLFilterDesigner:
|
|
|
|
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
|
|
|
|
|
2020-01-20 18:46:41 +00:00
|
|
|
def __init__(self, fs):
|
|
|
|
"""Initialize a filter bank designer.
|
|
|
|
|
|
|
|
Args:
|
|
|
|
fs: Sampling frequency [Hz]
|
|
|
|
"""
|
|
|
|
self.fs = fs
|
|
|
|
|
2020-01-17 21:10:42 +00:00
|
|
|
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 A(self, f):
|
|
|
|
"""
|
|
|
|
Computes the linear A-weighting freqency response. Hence, to obtain
|
|
|
|
A-weighted values, the *amplitude* need to be multiplied with this value.
|
|
|
|
Hence, to correct dB levels, the value of 20*log(A) needs to be added to
|
|
|
|
the level
|
|
|
|
|
|
|
|
Args:
|
|
|
|
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_uncor(self, f):
|
|
|
|
"""
|
|
|
|
Computes the uncorrected frequency response of the C-filter
|
|
|
|
"""
|
|
|
|
fsq = f**2
|
|
|
|
num = self.f4sq*fsq
|
|
|
|
denom1 = (fsq+self.f1**2)
|
|
|
|
denom2 = (fsq+self.f4**2)
|
|
|
|
return num/(denom1*denom2)
|
|
|
|
|
|
|
|
|
|
|
|
def C(self, f):
|
|
|
|
"""
|
|
|
|
Computes the linear A-weighting freqency response
|
|
|
|
"""
|
|
|
|
Cuncor = self._C_uncor(f)
|
|
|
|
C1000 = self._C_uncor(self.fr)
|
|
|
|
return Cuncor/C1000
|
|
|
|
|
|
|
|
|
2020-01-20 18:46:41 +00:00
|
|
|
def A_fir_design(self):
|
|
|
|
|
|
|
|
fs = self.fs
|
2020-01-17 21:10:42 +00:00
|
|
|
|
|
|
|
assert int(fs) == 48000
|
|
|
|
freq_design = np.linspace(0, 17e3, 3000)
|
|
|
|
freq_design[-1] = fs/2
|
|
|
|
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
|
|
|
|
|
|
|
|
|
2020-01-20 18:46:41 +00:00
|
|
|
def C_fir_design(self):
|
|
|
|
fs = self.fs
|
2020-01-17 21:10:42 +00:00
|
|
|
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
|
|
|
|
|
2020-01-20 18:46:41 +00:00
|
|
|
def C_Sos_design(self):
|
2020-01-17 21:10:42 +00:00
|
|
|
"""
|
|
|
|
Create filter coefficients of the C-weighting filter. Uses the bilinear
|
|
|
|
transform to convert the analog filter to a digital one.
|
|
|
|
|
|
|
|
Returns:
|
|
|
|
Sos: Second order sections
|
|
|
|
"""
|
2020-01-20 18:46:41 +00:00
|
|
|
fs = self.fs
|
2020-01-17 21:10:42 +00:00
|
|
|
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
|
2020-01-20 18:46:41 +00:00
|
|
|
|
|
|
|
def A_Sos_design(self):
|
|
|
|
"""
|
|
|
|
Create filter coefficients of the A-weighting filter. Uses the bilinear
|
|
|
|
transform to convert the analog filter to a digital one.
|
|
|
|
|
|
|
|
Returns:
|
|
|
|
Sos: Second order sections
|
|
|
|
"""
|
|
|
|
fs = self.fs
|
|
|
|
p1 = 2*np.pi*self.f1
|
|
|
|
p2 = 2*np.pi*self.f2
|
|
|
|
p3 = 2*np.pi*self.f3
|
|
|
|
p4 = 2*np.pi*self.f4
|
|
|
|
zeros_analog = [0,0,0,0]
|
|
|
|
poles_analog = [p1,p1, p2,p3,p4, p4]
|
|
|
|
k_analog = p4**2/self._A_uncor(self.fr)
|
|
|
|
|
|
|
|
z, p, k = bilinear_zpk(zeros_analog, poles_analog, k_analog, fs)
|
|
|
|
sos = zpk2sos(z, p, k)
|
|
|
|
return sos
|
|
|
|
|
2020-01-18 12:45:17 +00:00
|
|
|
|
|
|
|
def show_Afir():
|
|
|
|
from asceefig.plot import Figure
|
|
|
|
|
|
|
|
fs = 48000.
|
|
|
|
freq_design = np.linspace(0, 17e3, 3000)
|
|
|
|
freq_design[-1] = fs/2
|
|
|
|
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='hann'))
|
|
|
|
firs.append(A_fir_design())
|
|
|
|
# from scipy.signal import iirdesign
|
|
|
|
# b,a = iirdesign()
|
|
|
|
freq_check = np.logspace(0, np.log10(fs/2), 5000)
|
|
|
|
f = Figure()
|
|
|
|
|
|
|
|
f.semilogx(freq_check, 20*np.log10(A(freq_check)))
|
|
|
|
for fir in firs:
|
|
|
|
H = freqResponse(fs, freq_check, fir)
|
|
|
|
f.plot(freq_check, 20*np.log10(np.abs(H)))
|
|
|
|
|
|
|
|
f.fig.get_axes()[0].set_ylim(-75, 3)
|
|
|
|
|
|
|
|
|
|
|
|
def show_Cfir():
|
|
|
|
from asceefig.plot import Figure
|
|
|
|
|
|
|
|
fs = 48000.
|
|
|
|
freq_design = np.linspace(0, 17e3, 3000)
|
|
|
|
freq_design[-1] = fs/2
|
|
|
|
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='hann'))
|
|
|
|
firs.append(C_fir_design())
|
|
|
|
# from scipy.signal import iirdesign
|
|
|
|
# b,a = iirdesign()
|
|
|
|
freq_check = np.logspace(0, np.log10(fs/2), 5000)
|
|
|
|
f = Figure()
|
|
|
|
|
|
|
|
f.semilogx(freq_check, 20*np.log10(C(freq_check)))
|
|
|
|
for fir in firs:
|
|
|
|
H = freqResponse(fs, freq_check, fir)
|
|
|
|
f.plot(freq_check, 20*np.log10(np.abs(H)))
|
|
|
|
|
|
|
|
f.fig.get_axes()[0].set_ylim(-30, 1)
|
2020-01-17 21:10:42 +00:00
|
|
|
|