ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 5 Pull Requests Actions Releases 14 Activity

Merge branch 'splweighting_sos' into slm_c_impl

Browse Source
This commit is contained in:
Anne de Jong 2020-01-18 13:45:28 +01:00
commit b5088bef14
2 changed files with 128 additions and 121 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

247
lasp/filter/soundpressureweighting.py
View File

@ -7,151 +7,158 @@ 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',
'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
__all__ = ['SPLFilterDesigner']
def A_uncor(f):
"""
Computes the uncorrected frequency response of the A-filter
"""
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))
class SPLFilterDesigner:
fr = 1000.
fL = 10**1.5
fH = 10**3.9
def A(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
fLsq = fL**2
fHsq = fH**2
frsq = fr**2
fA = 10**2.45
D = np.sqrt(.5)
Args:
f: Frequency array to compute values for
Returns:
A(f) for each frequency
"""
Auncor = A_uncor(f)
A1000 = A_uncor(1000.)
return Auncor/A1000
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(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):
"""
Computes the uncorrected frequency response of the C-filter
"""
fsq = f**2
num = f4sq*fsq
denom1 = (fsq+f1**2)
denom2 = (fsq+f4**2)
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(f):
"""
Computes the linear A-weighting freqency response
"""
Cuncor = C_uncor(f)
C1000 = C_uncor(1000.)
return Cuncor/C1000
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 A_fir_design():
fs = 48000.
freq_design = np.linspace(0, 17e3, 3000)
freq_design[-1] = fs/2
amp_design = 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(self, f):
"""
Computes the linear A-weighting freqency response
"""
Cuncor = self._C_uncor(f)
C1000 = self._C_uncor(self.fr)
return Cuncor/C1000
def C_fir_design():
fs = 48000.
freq_design = np.linspace(0, 17e3, 3000)
freq_design[-1] = fs/2
amp_design = C(freq_design)
amp_design[-1] = 0.
def A_fir_design(self, fs):
L = 2048 # Filter order
fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
window='rectangular')
return fir
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
def show_Afir():
from asceefig.plot import Figure
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.
fs = 48000.
freq_design = np.linspace(0, 17e3, 3000)
freq_design[-1] = fs/2
amp_design = A(freq_design)
amp_design[-1] = 0.
firs = []
L = 2048 # Filter order
fir = arbitrary_fir_design(fs, L, freq_design, amp_design,
window='rectangular')
return fir
# 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()
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.
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)))
Args:
fs: Sampling frequency [Hz]
f.fig.get_axes()[0].set_ylim(-75, 3)
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
def show_Cfir():
from asceefig.plot import Figure
def A_Sos_design(self, fs):
"""
Create filter coefficients of the A-weighting filter. Uses the bilinear
transform to convert the analog filter to a digital one.
fs = 48000.
freq_design = np.linspace(0, 17e3, 3000)
freq_design[-1] = fs/2
amp_design = C(freq_design)
amp_design[-1] = 0.
firs = []
Args:
fs: Sampling frequency [Hz]
# 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()
Returns:
Sos: Second order sections
"""
# Poles of A-filter
p1 = 2*np.pi*self.f1
p2 = 2*np.pi*self.f2
p3 = 2*np.pi*self.f3
p4 = 2*np.pi*self.f4
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)))
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
f.fig.get_axes()[0].set_ylim(-30, 1)

2
lasp/lasp_weighcal.py
View File

@ -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