diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..66a6eb8 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -108,6 +108,7 @@ from . import utils from . import rooms from . import iso3741 +from . import iec61672 from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +131,7 @@ 'generate', 'utils', 'iso3741', + 'iec61672', 'rooms', # Functions diff --git a/pytta/iec61672.py b/pytta/iec61672.py new file mode 100644 index 0000000..fdce5fc --- /dev/null +++ b/pytta/iec61672.py @@ -0,0 +1,113 @@ +""" +PyTTa IEC 61672-1 sound level meter: +----------------------------------- +Frequency weightings A, C and Z, time weightings F and S, and the levels +of a sound level meter of IEC 61672-1:2013, plus the percentile levels of a +time-weighted level. + + >>> pytta.iec61672.leq(signalObj, weighting='A') # LAeq, dB + +Levels are in dB re the dBRef of each channel (20 uPa for a channel in Pa) +and are returned with one value per channel. +""" +import numpy as np +import scipy.signal as ss + +_TIME_CONSTANTS = {'F': 0.125, 'S': 1.0} # 5.8.1 + + +def _poles(): + """Pole frequencies f1 to f4 in Hz, from Annex E (E.2.4, E.3.4).""" + fr, fL, fH, D = 1000.0, 10 ** 1.5, 10 ** 3.9, 1 / np.sqrt(2) + b = (fr ** 2 + fL ** 2 * fH ** 2 / fr ** 2 - D * (fL ** 2 + fH ** 2)) \ + / (1 - D) + c = fL ** 2 * fH ** 2 + f1 = np.sqrt((-b - np.sqrt(b ** 2 - 4 * c)) / 2) + f4 = np.sqrt((-b + np.sqrt(b ** 2 - 4 * c)) / 2) + fA = 10 ** 2.45 + return f1, (3 - 5 ** 0.5) / 2 * fA, (3 + 5 ** 0.5) / 2 * fA, f4 + + +def weighting_filter(kind, fs): + """ + Second-order sections of the frequency weighting 'A', 'C' or 'Z' at + the sampling rate fs: the poles and zeros of Annex E through the + bilinear transform, with 0 dB at 1 kHz. 'Z' returns None (no filter). + """ + if kind not in ('A', 'C', 'Z'): + raise ValueError("weighting must be 'A', 'C' or 'Z'.") + if kind == 'Z': + return None + f1, f2, f3, f4 = _poles() + poles = [f1, f1, f4, f4] + ([f2, f3] if kind == 'A' else []) + zeros = [0.0] * (4 if kind == 'A' else 2) + sos = ss.zpk2sos(*ss.bilinear_zpk(zeros, [-2 * np.pi * p for p in poles], + 1.0, fs)) + _, h = ss.sosfreqz(sos, worN=np.array([1000.0]), fs=fs) + sos[0, :3] /= abs(h[0]) + return sos + + +def _weighted(signalObj, weighting): + """Frequency-weighted samples (samples, channels) and the dB reference + of each channel.""" + x = np.asarray(signalObj.timeSignal, dtype=float) + sos = weighting_filter(weighting, signalObj.samplingRate) + if sos is not None: + x = ss.sosfilt(sos, x, axis=0) + refs = np.array([signalObj.channels[ch].dBRef + for ch in signalObj.channels.mapping], dtype=float) + return x, refs + + +def time_weighted_level(signalObj, weighting='A', time='F'): + """ + Exponentially time-weighted sound level of 5.8 (F: 0.125 s, S: 1 s), + at every sample: (time in s, level in dB of shape (samples, channels)). + """ + if time not in _TIME_CONSTANTS: + raise ValueError("time must be 'F' or 'S'.") + x, refs = _weighted(signalObj, weighting) + alpha = np.exp(-1 / (_TIME_CONSTANTS[time] * signalObj.samplingRate)) + meanSquare = ss.lfilter([1 - alpha], [1, -alpha], x ** 2, axis=0) + with np.errstate(divide='ignore'): + level = 10 * np.log10(meanSquare / refs ** 2) + return np.arange(x.shape[0]) / signalObj.samplingRate, level + + +def leq(signalObj, weighting='A'): + """Time-averaged (equivalent continuous) sound level over the signal.""" + x, refs = _weighted(signalObj, weighting) + return 10 * np.log10(np.mean(x ** 2, axis=0) / refs ** 2) + + +def sel(signalObj, weighting='A'): + """Sound exposure level, the time integral of the squared weighted + signal re the reference squared times 1 s.""" + x, refs = _weighted(signalObj, weighting) + return 10 * np.log10(np.sum(x ** 2, axis=0) / signalObj.samplingRate + / refs ** 2) + + +def lmax(signalObj, weighting='A', time='F'): + """Maximum time-weighted sound level.""" + return np.max(time_weighted_level(signalObj, weighting, time)[1], axis=0) + + +def lpeak(signalObj, weighting='C'): + """Peak sound level, from the largest magnitude of the weighted + signal (C by default, as LCpeak).""" + x, refs = _weighted(signalObj, weighting) + return 20 * np.log10(np.max(np.abs(x), axis=0) / refs) + + +def percentile_levels(signalObj, percents=(10, 50, 90), weighting='A', + time='F'): + """ + Percentile levels L_N: the time-weighted level exceeded during N % of + the signal, for each N in percents. Returns an array of shape + (len(percents), channels). + """ + level = time_weighted_level(signalObj, weighting, time)[1] + return np.percentile(level, 100 - np.asarray(percents, dtype=float), + axis=0) diff --git a/tests/test_iec61672.py b/tests/test_iec61672.py new file mode 100644 index 0000000..c6ec118 --- /dev/null +++ b/tests/test_iec61672.py @@ -0,0 +1,112 @@ +""" +Tests for pytta.iec61672, the sound level meter measures of IEC 61672-1:2013 +(frequency weightings A, C and Z, time weightings F and S, Leq, SEL, Lmax, +Lpeak) and the percentile levels of a time-weighted level. Expected values +and acceptance limits are those of the standard. Each test pins one rule. +""" +import os +import unittest + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import numpy as np +import scipy.signal as ss + +import pytta +from pytta import iec61672 + +FS = 48000 +P94 = np.sqrt(2) # sine amplitude of 94 dB SPL, in Pa + +# Table 3: nominal frequency, A, C, class 1 acceptance limits (+, -) +TABLE3 = [(10, -70.4, -14.3, 3.0, np.inf), (16, -56.7, -8.5, 2.0, 4.0), + (31.5, -39.4, -3.0, 1.5, 1.5), (63, -26.2, -0.8, 1.0, 1.0), + (125, -16.1, -0.2, 1.0, 1.0), (250, -8.6, 0.0, 1.0, 1.0), + (500, -3.2, 0.0, 1.0, 1.0), (1000, 0.0, 0.0, 0.7, 0.7), + (2000, 1.2, -0.2, 1.0, 1.0), (4000, 1.0, -0.8, 1.0, 1.0), + (8000, -1.1, -3.0, 1.5, 2.5), (10000, -2.5, -4.4, 2.0, 3.0), + (12500, -4.3, -6.2, 2.0, 5.0), (16000, -6.6, -8.5, 2.5, 16.0)] + + +def in_pascal(signalObj): + for ch in signalObj.channels.mapping: + signalObj.channels[ch].unit = 'Pa' + return signalObj + + +def sine(seconds, amplitude, freq=1000.0): + t = np.arange(int(round(seconds * FS))) / FS + return amplitude * np.sin(2 * np.pi * freq * t) + + +class TestIEC61672(unittest.TestCase): + + def test_frequency_weightings_within_class_1(self): + # At the exact frequencies of Table 3 (base-ten, 1000 * 10^(n/10)) + # the response deviates from the design goals by less than the + # class 1 acceptance limits, at 48 kHz. + for kind, column in (('A', 1), ('C', 2)): + sos = iec61672.weighting_filter(kind, FS) + for row in TABLE3: + f = 1000 * 10 ** (round(10 * np.log10(row[0] / 1000)) / 10) + _, h = ss.sosfreqz(sos, worN=np.array([f]), fs=FS) + dev = 20 * np.log10(abs(h[0])) - row[column] + self.assertLessEqual(dev, row[3] + 0.05, (kind, row[0])) + self.assertGreaterEqual(dev, -row[4] - 0.05, (kind, row[0])) + + def test_toneburst_response_of_time_weighting_f(self): + # Table 4: LAFmax - LA of single 4 kHz tonebursts, and LAE - LA, + # which is 10 lg(Tb / 1 s). + table4 = [(1000, 0.0, 0.5, 0.5), (500, -0.1, 0.5, 0.5), + (200, -1.0, 0.5, 0.5), (100, -2.6, 1.0, 1.0), + (50, -4.8, 1.0, 1.0), (20, -8.3, 1.0, 1.0), + (10, -11.1, 1.0, 1.0), (5, -14.1, 1.0, 1.0), + (2, -18.0, 1.0, 1.5)] + steady = pytta.SignalObj(sine(2, P94, 4000), 'time', FS) + LA = iec61672.leq(in_pascal(steady))[0] + for ms, ref, up, down in table4: + burst = np.zeros(2 * FS) + burst[FS // 2:FS // 2 + int(ms * FS / 1000)] = \ + sine(ms / 1000, P94, 4000) + sig = in_pascal(pytta.SignalObj(burst, 'time', FS)) + response = iec61672.lmax(sig, time='F')[0] - LA + self.assertLessEqual(response, ref + up, ms) + self.assertGreaterEqual(response, ref - down, ms) + self.assertAlmostEqual(iec61672.sel(sig)[0] - LA, + 10 * np.log10(ms / 1000), delta=0.05) + + def test_levels_of_a_steady_94_db_tone(self): + # 1 Pa rms at 1 kHz: Leq and Lmax 94 dB with A (0 dB at 1 kHz), SEL + # of 10 s 104 dB, LCpeak of the amplitude sqrt(2) Pa 97.0 dB. + sig = in_pascal(pytta.SignalObj(sine(10, P94), 'time', FS)) + self.assertAlmostEqual(iec61672.leq(sig)[0], 94.0, delta=0.05) + self.assertAlmostEqual(iec61672.lmax(sig, time='S')[0], 94.0, + delta=0.05) + self.assertAlmostEqual(iec61672.sel(sig)[0], 104.0, delta=0.05) + # The peak is that of the steady tone: a 50 ms fade-in keeps the + # onset transient of the C filter out of it. + x = sine(2, P94) + x[:FS // 20] *= np.sin(np.linspace(0, np.pi / 2, FS // 20)) ** 2 + faded = in_pascal(pytta.SignalObj(x, 'time', FS)) + self.assertAlmostEqual(iec61672.lpeak(faded)[0], + 94 + 20 * np.log10(np.sqrt(2)), delta=0.05) + + def test_percentile_levels(self): + # 30 s at 60 dB then 30 s at 80 dB: L10 is 80 dB and L90 60 dB. + x = np.concatenate((sine(30, P94 / 10 ** (34 / 20)), + sine(30, P94 / 10 ** (14 / 20)))) + sig = in_pascal(pytta.SignalObj(x, 'time', FS)) + L10, L90 = iec61672.percentile_levels(sig, (10, 90))[:, 0] + self.assertAlmostEqual(L10, 80.0, delta=0.1) + self.assertAlmostEqual(L90, 60.0, delta=0.1) + + def test_rejects_unknown_weightings(self): + sig = pytta.SignalObj(np.zeros(FS), 'time', FS) + with self.assertRaises(ValueError): + iec61672.leq(sig, weighting='B') + with self.assertRaises(ValueError): + iec61672.lmax(sig, time='I') + + +if __name__ == '__main__': + unittest.main()