From 1df4465bcda63223ab24b3a455735bba3dc4d347 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 21:39:39 +0100 Subject: [PATCH] [FEATURE] add pytta.psychoacoustics with loudness, sharpness and roughness Three metrics of a one-channel SignalObj in Pa, each returning a dict of floats and NumPy arrays: - loudness_zwicker: ISO 532-1:2017, stationary and time-varying methods, free or diffuse field. N (sone), LN (phon), the specific loudness over 240 steps of 0.1 Bark and, time-varying, N(t), LN(t) and N'(t) every 2 ms with Nmax and N5. Port of SQAT's Loudness_ISO532_1 (BSD-3, Ella Manor's AARAE code after the C code of the standard). Signals not at 48 kHz are resampled as MATLAB's resample does (matches it to 2e-15). - sharpness_din45692: DIN 45692:2009, Formula (1) on that specific loudness, k = 0.11; S, or S(t) and S50 for the time-varying method. Written from the standard. - roughness_daniel1997: Daniel and Weber (1997), R(t), R'(t) over 47 half-Bark channels and Rmean (asper). Port of SQAT's Roughness_Daniel1997 (MIT, Dik Hermes et al.). The BSD-3 and MIT notices of the ported files are in the module header. Validation. Against SQAT on golden data from tests/golden/make_loudness_sqat.m and make_roughness_sqat.m (MATLAB R2026b, same float32 samples): stationary loudness within 3e-13, time- varying within 3.1e-8 sone (the nonlinear decay amplifies 3e-14 of filter rounding on the decay of a burst), roughness within 2.8e-5 asper (FFT rounding, 32 kHz input). ISO 532-1 Annex B: signal 1 (band levels) 83.2957 sone against 83.296; the standard's files of signals 2 to 5 within 0.034 sone (0.26 %); signals 6 to 25 keep N(t) within 5 % or 0.1 sone at every frame, max deviation 0.068 sone. DIN 45692: the 41 test signals of Tables A.2 and A.3 within 0.45 of the tolerance (worst 10.5 kHz band, -0.19 acum), reference noise 1.02 acum. 1 kHz at 40 dB gives 0.999976 sone; the reference AM tone 1.000001 asper. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/index.rst | 5 + docs/pyttadoc/psychoacoustics.rst | 7 + pytta/__init__.py | 6 + pytta/psychoacoustics.py | 967 +++++++++++++++++++++++++++++ tests/golden/loudness_sqat.mat | Bin 0 -> 40814 bytes tests/golden/make_loudness_sqat.m | 122 ++++ tests/golden/make_roughness_sqat.m | 82 +++ tests/golden/roughness_sqat.mat | Bin 0 -> 3725 bytes tests/test_psychoacoustics.py | 296 +++++++++ 9 files changed, 1485 insertions(+) create mode 100644 docs/pyttadoc/psychoacoustics.rst create mode 100644 pytta/psychoacoustics.py create mode 100644 tests/golden/loudness_sqat.mat create mode 100644 tests/golden/make_loudness_sqat.m create mode 100644 tests/golden/make_roughness_sqat.m create mode 100644 tests/golden/roughness_sqat.mat create mode 100644 tests/test_psychoacoustics.py diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..ab5f549 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -57,6 +57,10 @@ Modules: * pytta.iso3741: calculations according to the standard; + * pytta.psychoacoustics: + loudness (ISO 532-1), sharpness (DIN 45692) and roughness (Daniel + and Weber 1997) of a sound pressure signal; + .. toctree:: classes/index utils/index @@ -64,4 +68,5 @@ Modules: functions generate rooms + psychoacoustics properties diff --git a/docs/pyttadoc/psychoacoustics.rst b/docs/pyttadoc/psychoacoustics.rst new file mode 100644 index 0000000..c79aab0 --- /dev/null +++ b/docs/pyttadoc/psychoacoustics.rst @@ -0,0 +1,7 @@ +.. currentmodule:: pytta + +Psychoacoustics +=============== + +.. automodule:: pytta.psychoacoustics + :members: diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..a483f56 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -72,6 +72,10 @@ * pytta.iso3741: calculations according to the standard; + * pytta.psychoacoustics: + loudness (ISO 532-1), sharpness (DIN 45692) and roughness (Daniel + and Weber 1997) of a sound pressure signal; + You can find usage examples in the toolbox's examples folder. @@ -108,6 +112,7 @@ from . import utils from . import rooms from . import iso3741 +from . import psychoacoustics from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +135,7 @@ 'generate', 'utils', 'iso3741', + 'psychoacoustics', 'rooms', # Functions diff --git a/pytta/psychoacoustics.py b/pytta/psychoacoustics.py new file mode 100644 index 0000000..58c4e01 --- /dev/null +++ b/pytta/psychoacoustics.py @@ -0,0 +1,967 @@ +""" +PyTTa psychoacoustics: +---------------------- + +Psychoacoustic metrics of a sound pressure signal. + +Available functions: + + >>> pytta.psychoacoustics.loudness_zwicker(signal) + >>> pytta.psychoacoustics.sharpness_din45692(signal) + >>> pytta.psychoacoustics.roughness_daniel1997(signal) + +* loudness_zwicker: Zwicker loudness, ISO 532-1:2017, stationary and + time-varying sounds, free or diffuse field; +* sharpness_din45692: sharpness, DIN 45692:2009, from that loudness; +* roughness_daniel1997: roughness after Daniel and Weber (1997). + +Every function takes a one-channel SignalObj whose time signal is the sound +pressure in pascal, and returns a dict of floats and NumPy arrays (keys and +shapes in each docstring). For a signal stored in another unit, scale it +first (for example with SignalObj.calib_pressure). + +Origin of the code +------------------ + +loudness_zwicker and roughness_daniel1997 are ports to Python of the MATLAB +implementations in SQAT, the Sound Quality Analysis Toolbox +(https://github.com/ggrecow/SQAT), by Gil Felix Greco, Roberto Merino-Martinez, +Alejandro Osses and Michael J. B. Lotinga. The ported files carry their own +permissive licences (the rest of SQAT is GPL-3.0 and was not used): + +* Loudness_ISO532_1.m and Do_OB13_ISO532_1.m (BSD-3-Clause): MATLAB code by + Ella Manor for AARAE (2015) after the C code of ISO 532-1:2017, Annex A.4; + adapted and validated for SQAT by Gil Felix Greco (2023, 2025); stand-alone + filter bank by Alejandro Osses (2023); modified by Sergio Aguirre and + Gil Felix Greco (2026) to mirror the C reference code. The copyright notice + and the licence are reproduced below. +* Roughness_Daniel1997.m, private/Get_Hweight_roughness.m, + private/Get_gzi_roughness.m and the utilities Terhardt_filterbank.m, + Terhardt_filterbank_params.m, Get_Bark.m, calculate_a0.m, get_statistics.m + and get_exceeded_value.m (MIT): first MATLAB implementation of the Daniel + and Weber model by J. E. Schrader at TU Eindhoven, supervised by Dik + Hermes, then maintained by Dik Hermes; adapted for SQAT by Gil Felix Greco + and Alejandro Osses (2014-2025), parallelised by Mike Lotinga (2025), + restructured by Sergio Aguirre and Gil Felix Greco (2026) after the + revised implementation of Dik Hermes (2025). The copyright notices and the + licence are reproduced below. +* The N5 percentile of loudness_zwicker follows get_exceeded_value.m of the + same utilities (MIT). + +The port follows the MATLAB code operation by operation, so that it gives +SQAT's numbers; SQAT's resampling (MATLAB resample) is reproduced as well. + +sharpness_din45692 is written from the text of DIN 45692:2009 (Formula (1)). + +References +---------- + +ISO 532-1:2017. Acoustics - Methods for calculating loudness - Part 1: +Zwicker method. + +DIN 45692:2009. Messtechnische Simulation der Hoerempfindung Schaerfe +(Measurement technique for the simulation of the auditory sensation of +sharpness). + +Daniel, P., & Weber, R. (1997). Psychoacoustical roughness: implementation of +an optimized model. Acustica, 83, 113-123. + +Felix Greco, G., Merino-Martinez, R., Osses, A., & Langer, S. C. (2023). SQAT: +a MATLAB-based toolbox for quantitative sound quality analysis. INTER-NOISE +and NOISE-CON Congress and Conference Proceedings, InterNoise23, Chiba, +Japan. doi:10.3397/IN_2023_1075 + +Licences of the ported code +--------------------------- + +Loudness (Loudness_ISO532_1.m, Do_OB13_ISO532_1.m), BSD-3-Clause: + + Copyright (c) <2015>, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are + met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER + OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +Roughness (Roughness_Daniel1997.m and the files it uses) and the N5 +percentile (get_exceeded_value.m), MIT: + + copyright (c) 1999-2023 Dik Hermes + copyright (c) 2025 Gil Felix Greco (get_statistics.m, get_exceeded_value.m) + SQAT utilities Terhardt_filterbank_params.m, Get_Bark.m, calculate_a0.m: + the SQAT authors (MIT, as stated in their headers) + + Permission is hereby granted, free of charge, to any person obtaining a + copy of this software and associated documentation files (the + "Software"), to deal in the Software without restriction, including + without limitation the rights to use, copy, modify, merge, publish, + distribute, sublicense, and/or sell copies of the Software, and to permit + persons to whom the Software is furnished to do so, subject to the + following conditions: + + The above copyright notice and this permission notice shall be included + in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN + NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, + DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR + OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE + USE OR OTHER DEALINGS IN THE SOFTWARE. + +""" + +import warnings +from math import gcd + +import numpy as np +import scipy.signal as ss + +__all__ = ['loudness_zwicker', 'sharpness_din45692', 'roughness_daniel1997'] + +_FIELDS = ('free', 'diffuse') +_METHODS = ('stationary', 'time_varying') + + +# --------------------------------------------------------------------------- +# Common helpers +# --------------------------------------------------------------------------- + +def _mono(signal): + """Time signal (float64, 1-D) and sampling rate of a one-channel SignalObj.""" + x = np.asarray(signal.timeSignal, dtype=float) + if x.ndim == 2: + if x.shape[1] != 1: + raise ValueError('psychoacoustic metrics take one channel at a ' + f'time; the signal has {x.shape[1]} channels.') + x = x[:, 0] + return x, signal.samplingRate + + +def _resample(x, fs_out, fs_in): + """ + Resample x from fs_in to fs_out as MATLAB's resample(x, p, q) does. + + Polyphase FIR with a Kaiser window (beta 5) over 10 periods of the + slower rate on each side, filter delay removed, output length + ceil(len(x) * p / q). Reproduces MATLAB to rounding (2e-15 on a unit + variance noise). + """ + g = gcd(int(fs_out), int(fs_in)) + p, q = int(fs_out) // g, int(fs_in) // g + pqmax = max(p, q) + L = 2 * 10 * pqmax + 1 + m = np.arange(L) - (L - 1) / 2 + h = np.sinc(m / pqmax) * ss.windows.kaiser(L, 5.0) + h = p * h / np.sum(h) + Lhalf = (L - 1) / 2 + Lx = len(x) + nz = int(np.floor(q - np.mod(Lhalf, q))) + h = np.concatenate([np.zeros(nz), h]) + Lhalf += nz + delay = int(np.floor(np.ceil(Lhalf) / q)) + Ly = int(np.ceil(Lx * p / q)) + nz1 = 0 + while np.ceil(((Lx - 1) * p + len(h) + nz1) / q) - delay < Ly: + nz1 += 1 + h = np.concatenate([h, np.zeros(nz1)]) + y = ss.upfirdn(h, x, p, q) + return y[delay:delay + Ly] + + +def _exceeded(values, percent): + """Value exceeded during percent % of the time (SQAT get_exceeded_value).""" + idx = int(np.floor((100 - percent) / 100 * values.size)) + return np.sort(values)[max(idx, 1) - 1] + + +def _skip_index(time, time_skip): + """First frame of the statistics: the frame closest to time_skip.""" + return int(np.argmin(np.abs(time - time_skip))) + + +# --------------------------------------------------------------------------- +# Loudness, ISO 532-1:2017 (Zwicker) +# --------------------------------------------------------------------------- + +# One-third-octave filter bank at 48 kHz, ISO 532-1 Tables A.1 and A.2: three +# second-order sections per band, numerators _OB_B, denominators +# [1, -2, 1] - [0, a1, a2] (offsets per band and section, below), gain per band. +_OB_B = ((1., 2., 1.), (1., 0., -1.), (1., -2., 1.)) +_OB_AD = np.array([ # per band: (a1, a2) of sections 1, 2, 3 + [-6.70260e-004, 6.59453e-004, -3.75071e-004, 3.61926e-004, -3.06523e-004, 2.97634e-004], + [-8.47258e-004, 8.30131e-004, -4.76448e-004, 4.55616e-004, -3.88773e-004, 3.74685e-004], + [-1.07210e-003, 1.04496e-003, -6.06567e-004, 5.73553e-004, -4.94004e-004, 4.71677e-004], + [-1.35836e-003, 1.31535e-003, -7.74327e-004, 7.22007e-004, -6.29154e-004, 5.93771e-004], + [-1.72380e-003, 1.65564e-003, -9.91780e-004, 9.08866e-004, -8.03529e-004, 7.47455e-004], + [-2.19188e-003, 2.08388e-003, -1.27545e-003, 1.14406e-003, -1.02976e-003, 9.40900e-004], + [-2.79386e-003, 2.62274e-003, -1.64828e-003, 1.44006e-003, -1.32520e-003, 1.18438e-003], + [-3.57182e-003, 3.30071e-003, -2.14252e-003, 1.81258e-003, -1.71397e-003, 1.49082e-003], + [-4.58305e-003, 4.15355e-003, -2.80413e-003, 2.28135e-003, -2.23006e-003, 1.87646e-003], + [-5.90655e-003, 5.22622e-003, -3.69947e-003, 2.87118e-003, -2.92205e-003, 2.36178e-003], + [-7.65243e-003, 6.57493e-003, -4.92540e-003, 3.61318e-003, -3.86007e-003, 2.97240e-003], + [-1.00023e-002, 8.29610e-003, -6.63788e-003, 4.55999e-003, -5.15982e-003, 3.75306e-003], + [-1.31230e-002, 1.04220e-002, -9.02274e-003, 5.73132e-003, -6.94543e-003, 4.71734e-003], + [-1.73693e-002, 1.30947e-002, -1.24176e-002, 7.20526e-003, -9.46002e-003, 5.93145e-003], + [-2.31934e-002, 1.64308e-002, -1.73009e-002, 9.04761e-003, -1.30358e-002, 7.44926e-003], + [-3.13292e-002, 2.06370e-002, -2.44342e-002, 1.13731e-002, -1.82108e-002, 9.36778e-003], + [-4.28261e-002, 2.59325e-002, -3.49619e-002, 1.43046e-002, -2.57855e-002, 1.17912e-002], + [-5.91733e-002, 3.25054e-002, -5.06072e-002, 1.79513e-002, -3.69401e-002, 1.48094e-002], + [-8.26348e-002, 4.05894e-002, -7.40348e-002, 2.24476e-002, -5.34977e-002, 1.85371e-002], + [-1.17018e-001, 5.08116e-002, -1.09516e-001, 2.81387e-002, -7.85097e-002, 2.32872e-002], + [-1.67714e-001, 6.37872e-002, -1.63378e-001, 3.53729e-002, -1.16419e-001, 2.93723e-002], + [-2.42528e-001, 7.98576e-002, -2.45161e-001, 4.43370e-002, -1.73972e-001, 3.70015e-002], + [-3.53142e-001, 9.96330e-002, -3.69163e-001, 5.53535e-002, -2.61399e-001, 4.65428e-002], + [-5.16316e-001, 1.24177e-001, -5.55473e-001, 6.89403e-002, -3.93998e-001, 5.86715e-002], + [-7.56635e-001, 1.55023e-001, -8.34281e-001, 8.58123e-002, -5.94547e-001, 7.43960e-002], + [-1.10165e+000, 1.91713e-001, -1.23939e+000, 1.05243e-001, -8.91666e-001, 9.40354e-002], + [-1.58477e+000, 2.39049e-001, -1.80505e+000, 1.28794e-001, -1.32500e+000, 1.21333e-001], + [-2.50630e+000, 1.42308e-001, -2.19464e+000, 2.76470e-001, -1.90231e+000, 1.47304e-001], +]) +_OB_GAIN = np.array([ + 4.30764e-011, 8.59340e-011, 1.71424e-010, 3.41944e-010, + 6.82035e-010, 1.36026e-009, 2.71261e-009, 5.40870e-009, + 1.07826e-008, 2.14910e-008, 4.28228e-008, 8.54316e-008, + 1.70009e-007, 3.38215e-007, 6.71990e-007, 1.33531e-006, + 2.65172e-006, 5.25477e-006, 1.03780e-005, 2.04870e-005, + 4.05198e-005, 7.97914e-005, 1.56511e-004, 3.04954e-004, + 5.99157e-004, 1.16544e-003, 2.27488e-003, 3.91006e-003, +]) +_OB_FC = np.array([10 ** (((i - 1) - 16) / 10.) * 1000 for i in range(1, 29)]) + +# Table A.3: ranges (RAP) and level reductions (DLL) of the bands below 315 Hz +_RAP = np.array([45, 55, 65, 71, 80, 90, 100, 120], dtype=float) +_DLL = np.array([ + [-32, -24, -16, -10, -5, 0, -7, -3, 0, -2, 0], + [-29, -22, -15, -10, -4, 0, -7, -2, 0, -2, 0], + [-27, -19, -14, -9, -4, 0, -6, -2, 0, -2, 0], + [-25, -17, -12, -9, -3, 0, -5, -2, 0, -2, 0], + [-23, -16, -11, -7, -3, 0, -4, -1, 0, -1, 0], + [-20, -14, -10, -6, -3, 0, -4, -1, 0, -1, 0], + [-18, -12, -9, -6, -2, 0, -3, -1, 0, -1, 0], + [-15, -10, -8, -4, -2, 0, -3, -1, 0, -1, 0]], dtype=float) +# Table A.5 (free to diffuse field), A.6 (threshold in quiet), A.7 (critical +# band adaptation) and A.4 (a0, transmission of the middle ear) +_DDF = np.array([0, 0, 0.5, 0.9, 1.2, 1.6, 2.3, 2.8, 3, 2, 0, -1.4, -2, -1.9, + -1, 0.5, 3, 4, 4.3, 4]) +_LTQ = np.array([30, 18, 12, 8, 7, 6, 5, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3], + dtype=float) +_DCB = np.array([-0.25, -0.6, -0.8, -0.8, -0.5, 0, 0.5, 1.1, 1.5, 1.7, 1.8, 1.8, + 1.7, 1.6, 1.4, 1.2, 0.8, 0.5, 0, -0.5]) +_A0 = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -1.6, -3.2, -5.4, -5.6, -4, + -1.5, 2, 5, 12]) +# Table A.8: upper limits of the approximated critical bands (Bark); the +# 0.0001 is added as in the C code of ISO 532-1 +_ZUP = [z + 0.0001 for z in (.9, 1.8, 2.8, 3.5, 4.4, 5.4, 6.6, 7.9, 9.2, 10.6, + 12.3, 13.8, 15.2, 16.7, 18.1, 19.3, 20.6, 21.8, + 22.7, 23.6, 24)] +# Table A.9: ranges of specific loudness (RNS) and steepness of the upper +# slopes (USL, one row per range, one column per critical band index 1..>7). +# Row 13, column 5 is 0.24 as in SQAT; the corrected table of ISO 532-1:2017 +# gives 0.22 there. Kept at 0.24 so that the results equal SQAT's. +_RNS = [21.5, 18, 15.1, 11.5, 9, 6.1, 4.4, 3.1, 2.13, 1.36, 0.82, 0.42, 0.30, + 0.22, 0.15, 0.10, 0.035, 0] +_USL = [ + [13, 8.2, 6.3, 5.5, 5.5, 5.5, 5.5, 5.5], + [9, 7.5, 6, 5.1, 4.5, 4.5, 4.5, 4.5], + [7.8, 6.7, 5.6, 4.9, 4.4, 3.9, 3.9, 3.9], + [6.2, 5.4, 4.6, 4.0, 3.5, 3.2, 3.2, 3.2], + [4.5, 3.8, 3.6, 3.2, 2.9, 2.7, 2.7, 2.7], + [3.7, 3.0, 2.8, 2.35, 2.2, 2.2, 2.2, 2.2], + [2.9, 2.3, 2.1, 1.9, 1.8, 1.7, 1.7, 1.7], + [2.4, 1.7, 1.5, 1.35, 1.3, 1.3, 1.3, 1.3], + [1.95, 1.45, 1.3, 1.15, 1.1, 1.1, 1.1, 1.1], + [1.5, 1.2, 0.94, 0.86, 0.82, 0.82, 0.82, 0.82], + [0.72, 0.67, 0.64, 0.63, 0.62, 0.62, 0.62, 0.62], + [0.59, 0.53, 0.51, 0.50, 0.42, 0.42, 0.42, 0.42], + [0.40, 0.33, 0.26, 0.24, 0.24, 0.22, 0.22, 0.22], + [0.27, 0.21, 0.20, 0.18, 0.17, 0.17, 0.17, 0.17], + [0.16, 0.15, 0.14, 0.12, 0.11, 0.11, 0.11, 0.11], + [0.12, 0.11, 0.10, 0.08, 0.08, 0.08, 0.08, 0.08], + [0.09, 0.08, 0.07, 0.06, 0.06, 0.06, 0.06, 0.05], + [0.06, 0.05, 0.03, 0.02, 0.02, 0.02, 0.02, 0.02]] + +_I_REF = 4e-10 # (2e-5 Pa)^2 +_TINY = 1e-12 # added to the intensities before the logarithm +_FS = 48000 # the filter bank is defined at this rate +_SR_LEVEL = 2000 # rate of the band levels (time-varying method) +_SR_LOUDNESS = 500 # rate of the output (time-varying method) +_NL_ITER = 24 # inner iterations of the nonlinear decay per level sample +_BARK = np.arange(1, 241) / 10 + + +def _third_octave_bands(x): + """Signals of the 28 one-third-octave bands (25 Hz to 12.5 kHz) at 48 kHz.""" + out = np.empty((x.size, 28)) + for n in range(28): + y = x + for s in range(3): + a = (1., -2. - _OB_AD[n, 2 * s], 1. - _OB_AD[n, 2 * s + 1]) + y = ss.lfilter(_OB_B[s], a, y) + out[:, n] = _OB_GAIN[n] * y + return out + + +def _band_levels(x, time_varying, time_skip): + """ + One-third-octave band levels (dB): one row for the stationary method + (mean square after time_skip), one row per 0.5 ms for the time-varying + method (squared, smoothed by three first-order low-passes, decimated). + """ + bands = _third_octave_bands(x) ** 2 + n = x.size + if not time_varying: + skip = int(np.floor(time_skip * _FS)) + if skip >= n: + raise ValueError(f'time_skip ({time_skip} s) is not shorter than ' + f'the signal ({n / _FS} s).') + ms = np.sum(bands[skip:n], axis=0) / (n - skip) + return 10 * np.log10((ms + _TINY) / _I_REF)[np.newaxis, :] + dec = _FS // _SR_LEVEL + levels = np.empty((int(np.ceil(n / dec)), 28)) + for i in range(28): + tau = 2 / (3 * _OB_FC[i]) if _OB_FC[i] <= 1000 else 2 / (3 * 1000.) + a1 = np.exp(-1 / (_FS * tau)) + band = bands[:, i] + for _ in range(3): + band = ss.lfilter([1 - a1], [1, -a1], band) + levels[:, i] = 10 * np.log10((band[::dec] + _TINY) / _I_REF) + return levels + + +def _core_loudness(levels, diffuse): + """Core loudness of the 21 critical bands (ISO 532-1, steps 4 to 7).""" + nt = levels.shape[0] + # Step 4: weighting of the 11 bands below 315 Hz (Table A.3) + intens = np.empty((nt, 11)) + for i in range(11): + lev = levels[:, i] + k = np.zeros(nt, dtype=int) + active = np.ones(nt, dtype=bool) + for kk in range(7): + active &= lev > _RAP[kk] - _DLL[kk, i] + k[active] += 1 + intens[:, i] = 10 ** ((lev + _DLL[k, i]) / 10) + # Step 5: levels of the first three critical bands + cbi = np.stack([np.sum(intens[:, 0:6], axis=1), np.sum(intens[:, 6:9], axis=1), + np.sum(intens[:, 9:11], axis=1)], axis=1) + lcb = np.zeros_like(cbi) + pos = cbi > 0 + lcb[pos] = 10 * np.log10(cbi[pos]) + # Step 6: core loudness of each critical band + le = levels[:, 8:28].copy() + le[:, 0:3] = lcb + le = le - _A0 + if diffuse: + le = le + _DDF + core = np.zeros((nt, 21)) + above = le > _LTQ + mp1 = 0.0635 * 10 ** (0.025 * _LTQ) + mp2 = ((1 - 0.25) + 0.25 * 10 ** (0.1 * (le - _DCB - _LTQ))) ** 0.25 - 1 + c = mp1 * mp2 + c[c <= 0] = 0 + core[:, :20][above] = c[above] + # Step 7: correction of the lowest critical band + corr = 0.4 + 0.32 * core[:, 0] ** 0.2 + corr[corr > 1] = 1 + core[:, 0] = core[:, 0] * corr + return core + + +def _nonlinear_decay(core): + """Step 8: nonlinear temporal decay of the core loudness (in place).""" + t_short, t_long, t_var = 0.005, 0.015, 0.075 + dt = 1 / (_SR_LEVEL * _NL_ITER) + P = (t_var + t_long) / (t_var * t_short) + Q = 1 / (t_short * t_var) + l1 = -P / 2 + np.sqrt(P * P / 4 - Q) + l2 = -P / 2 - np.sqrt(P * P / 4 - Q) + den = t_var * (l1 - l2) + e1 = np.exp(l1 * dt) + e2 = np.exp(l2 * dt) + b1 = (e1 - e2) / den + b2 = ((t_var * l2 + 1) * e1 - (t_var * l1 + 1) * e2) / den + b3 = ((t_var * l1 + 1) * e1 - (t_var * l2 + 1) * e2) / den + b4 = (t_var * l1 + 1) * (t_var * l2 + 1) * (e1 - e2) / den + b5 = np.exp(-dt / t_long) + b6 = np.exp(-dt / t_var) + nt = core.shape[0] + for i in range(21): + col = core[:, i].tolist() + uo_last = 0.0 + u2_last = 0.0 + for j in range(nt): + ui = col[j] + if j < nt - 1: + delta = (col[j + 1] - ui) / _NL_ITER + n_inner = _NL_ITER + else: + delta = 0.0 + n_inner = 1 + for k in range(n_inner): + if ui < uo_last: + if uo_last > u2_last: + u2 = uo_last * b1 - u2_last * b2 + uo = uo_last * b3 - u2_last * b4 + if uo < ui: + uo = ui + if u2 > uo: + u2 = uo + else: + uo = uo_last * b5 + if uo < ui: + uo = ui + u2 = uo + elif abs(ui - uo_last) < 1e-5: + uo = ui + if uo > u2_last: + u2 = (u2_last - ui) * b6 + ui + else: + u2 = ui + else: + uo = ui + u2 = (u2_last - ui) * b6 + ui + uo_last = uo + u2_last = u2 + if k == 0: + col[j] = uo + ui += delta + core[:, i] = col + return core + + +def _slopes(core): + """ + Step 9: total loudness and specific loudness (240 values, 0.1 Bark) of + one row of core loudness, with the upper slopes of Table A.9. + """ + N = 0.0 + z1 = 0.0 + n1 = 0.0 + iz = 0 + z = 0.1 + j = 17 + ns = [0.0] * 245 # 240 used; slack for the float accumulation of k + for i in range(21): + ig = min(i, 8) - 1 # column of _USL; unused for i = 0 + zup = _ZUP[i] + cl = core[i] + while z1 < zup: + if n1 <= cl: + if n1 < cl: + j = 0 + while _RNS[j] >= cl and j < 17: + j += 1 + z2 = zup + n2 = cl + N += n2 * (z2 - z1) + k = z + while k <= z2: + ns[iz] = n2 + iz += 1 + k += 1 / 10 + z = k + else: + n2 = _RNS[j] + if n2 < cl: + n2 = cl + usl = _USL[j][ig] + dz = (n1 - n2) / usl + z2 = z1 + dz + if z2 > zup: + z2 = zup + dz = z2 - z1 + n2 = n1 - dz * usl + N += dz * (n1 + n2) / 2 + k = z + while k <= z2: + ns[iz] = n1 - (k - z1) * usl + iz += 1 + k += 1 / 10 + z = k + while n2 <= _RNS[j] and j < 17: + j += 1 + z1 = z2 + n1 = n2 + return max(N, 0.0), ns[:240] + + +def _loudness_level(N): + """Loudness level (phon) of a loudness N (sone), as in ISO 532-1.""" + N = np.asarray(N, dtype=float) + LN = 40 * N ** 0.35 + big = N >= 1 + LN[big] = 40 + 10 * np.log2(N[big]) + LN[LN < 3] = 3 + return LN + + +def _temporal_weighting(N): + """Step 10: two first-order low-passes (3.5 ms and 70 ms) combined.""" + dec = _FS // _SR_LEVEL + N = np.asarray(N, dtype=float).tolist() + out = [] + for tau in (3.5e-3, 70e-3): + a1 = np.exp(-1 / (_SR_LEVEL * dec * tau)) + b0 = 1 - a1 + y1 = 0.0 + y = [] + for i in range(len(N)): + x0 = N[i] + y1 = b0 * x0 + a1 * y1 + y.append(y1) + if i < len(N) - 1: + xd = (N[i + 1] - x0) / dec + for _ in range(dec - 1): + x0 = x0 + xd + y1 = b0 * x0 + a1 * y1 + out.append(np.array(y)) + return 0.47 * out[0] + 0.53 * out[1] + + +def _loudness_from_third_octave_levels(levels, field='free'): + """ + Stationary loudness from the 28 one-third-octave band levels (dB SPL, + 25 Hz to 12.5 kHz), ISO 532-1, 5.2 (the method 0 of SQAT). + + Returns a dict as loudness_zwicker with method='stationary'. + """ + levels = np.asarray(levels, dtype=float).reshape(1, 28) + core = _core_loudness(levels, field == 'diffuse') + N, ns = _slopes(core[0].tolist()) + return {'N': N, 'LN': float(_loudness_level([N])[0]), + 'specific': np.array(ns), 'bark': _BARK.copy()} + + +def loudness_zwicker(signal, field='free', method='stationary', time_skip=0.0): + """ + Loudness after Zwicker, ISO 532-1:2017. + + Port of Loudness_ISO532_1 of SQAT (see the module docstring); results + equal SQAT's to rounding. A signal at a rate other than 48 kHz is + resampled to 48 kHz first, as in SQAT (MATLAB resample). + + Parameters (default), (type): + ----------------------------- + + * signal (), (SignalObj): + One channel, sound pressure in Pa; + + * field ('free'), (str): + 'free' or 'diffuse' sound field; + + * method ('stationary'), (str): + 'stationary' (ISO 532-1, clause 5) or 'time_varying' (clause 6); + + * time_skip (0.0), (float): + Seconds skipped at the start. Stationary method: excluded from + the band levels (must be shorter than the signal). Time-varying + method: excluded from Nmax and N5 only. + + Returns: + -------- + + dict. Stationary method: + + * 'N': loudness, sone (float); + * 'LN': loudness level, phon (float); + * 'specific': specific loudness N', sone/Bark, shape (240,); + * 'bark': critical band rate of 'specific', 0.1 to 24 Bark in + steps of 0.1, shape (240,). + + Time-varying method, nt output frames every 2 ms: + + * 'N': loudness N(t), sone, shape (nt,); + * 'LN': loudness level LN(t), phon, shape (nt,); + * 'specific': specific loudness N'(t), sone/Bark, shape (nt, 240); + * 'time': time of each frame, s, shape (nt,); + * 'bark': as above, shape (240,); + * 'Nmax', 'N5': maximum loudness and loudness exceeded 5 % of + the time from time_skip on, sone (float). + + """ + if field not in _FIELDS: + raise ValueError(f"field must be one of {_FIELDS}, got {field!r}.") + if method not in _METHODS: + raise ValueError(f"method must be one of {_METHODS}, got {method!r}.") + x, fs = _mono(signal) + if fs != _FS: + x = _resample(x, _FS, fs) + time_varying = method == 'time_varying' + levels = _band_levels(x, time_varying, time_skip) + core = _core_loudness(levels, field == 'diffuse') + if not time_varying: + N, ns = _slopes(core[0].tolist()) + return {'N': N, 'LN': float(_loudness_level([N])[0]), + 'specific': np.array(ns), 'bark': _BARK.copy()} + + core = _nonlinear_decay(core) + dec = _SR_LEVEL // _SR_LOUDNESS + rows = core.tolist() + N_level = np.empty(len(rows)) + spec = [] + for l, row in enumerate(rows): + N_level[l], ns = _slopes(row) + if l % dec == 0: + spec.append(ns) + N_t = _temporal_weighting(N_level)[::dec] + time = np.arange(N_t.size) * 2e-3 + tail = N_t[_skip_index(time, time_skip):] + return {'N': N_t, 'LN': _loudness_level(N_t), 'specific': np.array(spec), + 'time': time, 'bark': _BARK.copy(), + 'Nmax': float(np.max(tail)), 'N5': float(_exceeded(tail, 5))} + + +# --------------------------------------------------------------------------- +# Sharpness, DIN 45692:2009 +# --------------------------------------------------------------------------- + +# Normalisation constant k of Formula (1), 0.105 <= k < 0.115, meant to bring +# the reference sound (narrowband noise 920 Hz to 1080 Hz, 60 dB) to 1.00 +# acum. k = 0.11 gives 1.02 to 1.03 acum for that noise synthesised with a flat +# spectrum (within the 0.05 acum tolerance of DIN 45692, 6) and keeps the 41 +# test signals of Tables A.2 and A.3 within 0.45 of their tolerance. The k +# that brings the synthetic noise to exactly 1.00 (0.1069) moves the +# broadband signals 0.1 acum down and the 10.5 kHz band out of tolerance. +_K_SHARPNESS = 0.11 + + +def _sharpness(specific): + """ + DIN 45692 Formula (1) over the 0.1 Bark grid of ISO 532-1 (the dz of + both integrals cancels). specific: (240,) or (nt, 240). Returns 0 where + the loudness is zero. + """ + z = _BARK + g = np.where(z <= 15.8, 1.0, 0.15 * np.exp(0.42 * (z - 15.8)) + 0.85) + num = np.sum(specific * g * z, axis=-1) + den = np.sum(specific, axis=-1) + with np.errstate(invalid='ignore', divide='ignore'): + S = np.where(den > 0, _K_SHARPNESS * num / den, 0.0) + return S + + +def sharpness_din45692(signal, field='free', method='stationary', time_skip=0.0): + """ + Sharpness, DIN 45692:2009. + + The first moment of the specific loudness of loudness_zwicker, weighted + by g(z) (1 up to 15.8 Bark, rising above), times k = 0.11 (Formula (1)). + The standard is written for the free field; a diffuse field result + should be reported as such. + + Parameters (default), (type): + ----------------------------- + + * signal (), (SignalObj): + One channel, sound pressure in Pa; + + * field ('free'), (str): + 'free' or 'diffuse', passed to loudness_zwicker; + + * method ('stationary'), (str): + 'stationary' or 'time_varying', passed to loudness_zwicker; + + * time_skip (0.0), (float): + Passed to loudness_zwicker; for the time-varying method it also + starts the frames of S50. + + Returns: + -------- + + dict. Stationary method: + + * 'S': sharpness, acum (float). + + Time-varying method: + + * 'S': sharpness S(t), acum, shape (nt,), one value per 2 ms + frame of the specific loudness (0 where the loudness is 0); + * 'time': time of each frame, s, shape (nt,); + * 'S50': median of S(t) from time_skip on, acum (float), the + single value DIN 45692 (clause 6) recommends for checks. + + """ + out = loudness_zwicker(signal, field=field, method=method, time_skip=time_skip) + S = _sharpness(out['specific']) + if method == 'stationary': + return {'S': float(S)} + tail = S[_skip_index(out['time'], time_skip):] + return {'S': S, 'time': out['time'], 'S50': float(np.median(tail))} + + +# --------------------------------------------------------------------------- +# Roughness, Daniel and Weber (1997) +# --------------------------------------------------------------------------- + +_CHNO = 47 # half-Bark channels +_R_WINDOW = 0.2 # s, analysis window (50 % overlap) +_R_DZ = 0.5 # Bark, integration step of the specific roughness +# Calibration: Daniel and Weber (1997), Eq. 9, cal = 0.25, re-derived in SQAT +# as 0.25 / R_ref with R_ref measured on SQAT's reference signal +_R_REF = 1.006197 + +# Critical bands (Get_Bark.m): band, lower edge (Hz), centre (Hz), centre (Bark) +_BARK_RAW = np.array([ + [0, 0, 50, 0.5], [1, 100, 150, 1.5], [2, 200, 250, 2.5], [3, 300, 350, 3.5], + [4, 400, 450, 4.5], [5, 510, 570, 5.5], [6, 630, 700, 6.5], + [7, 770, 840, 7.5], [8, 920, 1000, 8.5], [9, 1080, 1170, 9.5], + [10, 1270, 1370, 10.5], [11, 1480, 1600, 11.5], [12, 1720, 1850, 12.5], + [13, 2000, 2150, 13.5], [14, 2320, 2500, 14.5], [15, 2700, 2900, 15.5], + [16, 3150, 3400, 16.5], [17, 3700, 4000, 17.5], [18, 4400, 4800, 18.5], + [19, 5300, 5800, 19.5], [20, 6400, 7000, 20.5], [21, 7700, 8500, 21.5], + [22, 9500, 10500, 22.5], [23, 12000, 13500, 23.5], + [24, 15500, 20000, 24.5]], dtype=float) +# Hearing threshold of the excitation (Bark, dB), Terhardt_filterbank_params.m +_HTRES = np.array([ + [0, 130], [0.01, 70], [0.17, 60], [0.8, 30], [1, 25], [1.5, 20], [2, 15], + [3.3, 10], [4, 8.1], [5, 6.3], [6, 5], [8, 3.5], [10, 2.5], [12, 1.7], + [13.3, 0], [15, -2.5], [16, -4], [17, -3.7], [18, -1.5], [19, 1.4], + [20, 3.8], [21, 5], [22, 7.5], [23, 15], [24, 48], [24.5, 60], [25, 130]]) +# a0 transmission factor of the outer and middle ear, free field (Bark, dB), +# Fastl and Zwicker (2007), Fig. 8.18 (calculate_a0.m, 'fastl2007') +_A0TAB = np.array([ + [0, 0], [10, 0], [12, 1.15], [13, 2.31], [14, 3.85], [15, 5.62], + [16, 6.92], [16.5, 7.38], [17, 6.92], [18, 4.23], [18.5, 2.31], [19, 0], + [20, -1.43], [21, -2.59], [21.5, -3.57], [22, -5.19], [22.5, -7.41], + [23, -11.3], [23.5, -20], [24, -40], [25, -130], [26, -999]]) +# g(z) of the specific roughness, table of Dik Hermes (2025), Get_gzi_roughness.m +_GR = np.array([0.15, 0.26, 0.38, 0.47, 0.54, 0.65, 0.76, 0.83, 0.90, 0.98, + 0.98, 0.90, 0.80, 0.70, 0.62, 0.54, 0.49, 0.43, 0.39, 0.35, + 0.30, 0.30, 0.30, 0.30, 0.30]) +# Modulation weighting prototypes H(fmod) (Hz, weight), Get_Hweight_roughness.m +_H2 = np.array([[0, 0], [17, 0.8], [23, 0.95], [25, 0.975], [32, 1], + [37, 0.975], [48, 0.9], [67, 0.8], [90, 0.7], [114, 0.6], + [171, 0.4], [206, 0.3], [247, 0.2], [294, 0.1], [358, 0]]) +_H5 = np.array([[0, 0], [32, 0.8], [43, 0.95], [56, 1], [69, 0.975], + [92, 0.9], [120, 0.8], [142, 0.7], [165, 0.6], [231, 0.4], + [277, 0.3], [331, 0.2], [397, 0.1], [502, 0]]) +_H16 = np.array([[0, 0], [23.5, 0.4], [34, 0.6], [47, 0.8], [56, 0.9], + [63, 0.95], [79, 1], [100, 0.975], [115, 0.95], [135, 0.9], + [159, 0.85], [172, 0.8], [194, 0.7], [215, 0.6], [244, 0.5], + [290, 0.4], [348, 0.3], [415, 0.2], [500, 0.1], [645, 0]]) +_H21 = np.array([[0, 0], [19, 0.4], [44, 0.8], [52.5, 0.9], [58, 0.95], + [75, 1], [101.5, 0.95], [114.5, 0.9], [132.5, 0.85], + [143.5, 0.8], [165.5, 0.7], [197.5, 0.6], [241, 0.5], + [290, 0.4], [348, 0.3], [415, 0.2], [500, 0.1], [645, 0]]) +_H42 = np.array([[0, 0], [15, 0.4], [41, 0.8], [49, 0.9], [53, 0.965], + [64, 0.99], [71, 1], [88, 0.95], [94, 0.9], [106, 0.85], + [115, 0.8], [137, 0.7], [180, 0.6], [238, 0.5], [290, 0.4], + [348, 0.3], [415, 0.2], [500, 0.1], [645, 0]]) + + +def _mround(v): + """MATLAB round (half away from zero) for non-negative values.""" + return int(np.floor(v + 0.5)) + + +def _roughness_params(N, fs): + """Bins, Bark scale, thresholds, a0, H(fmod) and g(z) for windows of N samples.""" + df = fs / N + n0 = _mround(20 / df) + 1 # 1-based bin of 20 Hz + ntop = _mround(20e3 / df) + 1 # 1-based bin of 20 kHz + n01 = n0 - 1 + bins = np.arange(n0, ntop + 1) - 1 # 0-based bins from 20 Hz to 20 kHz + freqs = bins * df + f_sorted = np.sort(np.concatenate([_BARK_RAW[:, 1], _BARK_RAW[:, 2]])) + z_sorted = np.sort(np.concatenate([_BARK_RAW[:, 0], _BARK_RAW[:, 3]])) + bark = np.zeros(N) # Bark of each 0-based bin + bark[bins] = np.interp(freqs, f_sorted, z_sorted) + min_exc = np.interp(bark[bins], _HTRES[:, 0], _HTRES[:, 1]) + cf = np.array([_mround(f / df) for f in _BARK_RAW[1:25, 1]]) - n01 + 1 + bf = np.array([_mround(f / df) for f in _BARK_RAW[0:25, 2]]) - n01 + 1 + min_bf = min_exc[np.sort(np.concatenate([bf, cf])) - 1] + a0 = np.ones(N) + a0[bins] = 10 ** (np.interp(bark[bins], _A0TAB[:, 0], _A0TAB[:, 1]) / 20) + hw = np.zeros((_CHNO, N)) + for proto, fmax, rows in ((_H2, 358, (2,)), (_H5, 502, (5,)), + (_H16, 645, (16,)), (_H21, 645, (21,)), + (_H42, 645, (42,))): + last = int(np.floor((fmax / fs) * N)) + k = np.arange(3, last + 1) # 1-based bins, DC bins 1 and 2 left out + for r in rows: + hw[r - 1, k - 1] = np.interp((k - 1) * fs / N, proto[:, 0], proto[:, 1]) + for lo, hi, src in ((1, 4, 2), (6, 15, 5), (17, 20, 16), (22, 41, 21), + (43, 47, 42)): + for r in range(lo, hi + 1): + if r != src: + hw[r - 1] = hw[src - 1] + gzi = np.sqrt(np.interp(np.arange(1, _CHNO + 1) / 2, np.arange(25), _GR)) + return {'N': N, 'n01': n01, 'bins': bins, 'freqs': freqs, 'bark': bark, + 'min_exc': min_exc, 'min_bf': min_bf, 'a0': a0, 'hw': hw, + 'gzi': gzi} + + +def _terhardt_excitation(spec, prm): + """ + Excitation of the 47 half-Bark channels in the time domain (Terhardt + 1979 slopes) from the calibrated spectrum of one window. Returns the + (47, N) excitation and the number of components whose upper slope was + clamped to zero. + """ + N = prm['N'] + lg = np.abs(spec[prm['bins']]) + with np.errstate(divide='ignore'): + ldb = 20 * np.log10(lg) + which = np.nonzero(ldb > prm['min_exc'])[0] + if which.size == 0: + return np.zeros((_CHNO, N)), 0 + lw = ldb[which] + zw = prm['bark'][which + prm['n01']] + steep = -24 - (230 / prm['freqs'][which]) + (0.2 * lw) + s2 = np.where(steep < 0, steep, 0.0) + npos = prm['min_bf'].size + z1 = np.floor(2 * zw).astype(int) + z2 = np.ceil(2 * zw).astype(int) + cols = np.arange(1, npos + 1) + m1 = cols[np.newaxis, :] <= z1[:, np.newaxis] + m2 = (cols[np.newaxis, :] >= z2[:, np.newaxis]) & (cols <= _CHNO) + zk = 0.5 * cols[np.newaxis, :] + zi = zw[:, np.newaxis] + li = lw[:, np.newaxis] + stemp = np.zeros((which.size, npos)) + stemp = np.where(m1, -27 * (zi - zk) + li, stemp) + stemp = np.where(m2, s2[:, np.newaxis] * (zk - zi) + li, stemp) + reach = (m1 | m2) & (stemp > prm['min_bf'][np.newaxis, :]) + slopes = np.zeros_like(stemp) + slopes[reach] = 10 ** (stemp[reach] / 20) + amp = np.zeros((which.size, _CHNO)) + lgw = lg[which] + for i in range(1, _CHNO + 1): + if i != 1: + amp[:, i - 1] = slopes[:, i - 2] / lgw + sel = z2 > i + amp[sel, i - 1] = slopes[sel, i] / lgw[sel] + amp[z2 == i, i - 1] = 1 + amp[z1 == i, i - 1] = 1 + etmp = np.zeros((_CHNO, N), dtype=complex) + etmp[:, which + prm['n01']] = amp.T * spec[which + prm['n01']] + return 2 * N * np.real(np.fft.ifft(etmp, axis=1)), int(np.sum(steep >= 0)) + + +def _specific_roughness(ei, prm): + """Specific roughness (asper/Bark) of the 47 channels of one window.""" + N = prm['N'] + env = np.abs(ei) + h0 = np.mean(env, axis=1) + hbp = 2 * np.real(np.fft.ifft(np.fft.fft(env - h0[:, np.newaxis], axis=1) + * prm['hw'], axis=1)) + rms = np.sqrt(np.mean(hbp ** 2, axis=1)) + mdept = np.zeros(_CHNO) + pos = h0 > 0 + mdept[pos] = np.minimum(rms[pos] / h0[pos], 1) + ki = np.zeros(_CHNO - 2) + for k in range(_CHNO - 2): + a = hbp[k] - np.sum(hbp[k]) / N + b = hbp[k + 2] - np.sum(hbp[k + 2]) / N + c12 = np.dot(a, b) / (N - 1) + den = np.sqrt((np.dot(b, b) / (N - 1)) * (np.dot(a, a) / (N - 1))) + if den > 0: + ki[k] = c12 / den + g = prm['gzi'] + ri = np.zeros(_CHNO) + ri[0] = (g[0] * mdept[0] * ki[0]) ** 2 + ri[1] = (g[1] * mdept[1] * ki[1]) ** 2 + for k in range(2, _CHNO - 2): + ri[k] = (g[k] * mdept[k] * ki[k - 2] * ki[k]) ** 2 + ri[_CHNO - 2] = (g[_CHNO - 2] * mdept[_CHNO - 2] * ki[_CHNO - 4]) ** 2 + ri[_CHNO - 1] = (g[_CHNO - 1] * mdept[_CHNO - 1] * ki[_CHNO - 3]) ** 2 + return (0.25 / _R_REF / _R_DZ) * ri + + +def roughness_daniel1997(signal, time_skip=0.0): + """ + Roughness after Daniel and Weber (1997). + + Port of Roughness_Daniel1997 of SQAT (see the module docstring); results + equal SQAT's within the rounding of the FFT. Blackman windows of 200 ms + with 50 % overlap; per window, the spectrum is weighted by the a0 + transmission factor, spread over 47 half-Bark channels with Terhardt's + slopes, and the modulation depth of each channel envelope (weighted by + H(fmod)) and the correlation of channels 1 Bark apart give the specific + roughness. 1 asper: 1 kHz tone at 60 dB, 100 % amplitude modulated at + 70 Hz. Signals at rates other than 40960, 44100 and 48000 Hz are + resampled to 48 kHz first, as in SQAT (MATLAB resample). + + Parameters (default), (type): + ----------------------------- + + * signal (), (SignalObj): + One channel, sound pressure in Pa, at least 300 ms long (one + window and one hop); + + * time_skip (0.0), (float): + Seconds excluded from Rmean at the start. + + Returns: + -------- + + dict, nt windows: + + * 'R': roughness R(t), asper, shape (nt,); + * 'specific': specific roughness R'(t), asper/Bark, shape (nt, 47); + * 'time': start of each window, s, shape (nt,) (SQAT's + convention: (first sample index, 1-based) / fs); + * 'bark': critical band rate of 'specific', 0.5 to 23.5 Bark, + shape (47,); + * 'Rmean': mean of R(t) from the window closest to time_skip + on, asper (float). + + """ + x, fs = _mono(signal) + if fs not in (44100, 40960, 48000): + x = _resample(x, 48000, fs) + fs = 48000 + N = _mround(fs * _R_WINDOW) + hop = N // 2 + nwin = int(np.floor((x.size - N) / hop)) + if nwin < 1: + raise ValueError(f'the signal ({x.size / fs:.3f} s) is shorter than ' + f'one window and one hop ({1.5 * _R_WINDOW} s).') + prm = _roughness_params(N, fs) + window = ss.windows.blackman(N, sym=False) + l_cal = 20 * np.log10(1 / 20e-6) - 20 * np.log10(np.sqrt(2)) + amp_cal = 10 ** (l_cal / 20) * 2 / (N * np.mean(window)) + R = np.zeros(nwin) + spec = np.zeros((nwin, _CHNO)) + time = np.zeros(nwin) + clamped = 0 + for w in range(nwin): + start = w * hop + frame = np.fft.fft(x[start:start + N] * window * amp_cal) * prm['a0'] + ei, nclamp = _terhardt_excitation(frame, prm) + clamped += nclamp > 0 + ri = _specific_roughness(ei, prm) + spec[w] = ri + R[w] = _R_DZ * np.sum(ri) + time[w] = (start + 1) / fs + if clamped: + warnings.warn(f'Terhardt upper slope clamped to zero in {clamped} of ' + f'{nwin} window(s): a component exceeds 120 + 1150/f dB. ' + 'The roughness there is an extrapolation; check that the ' + 'signal is in Pa.', stacklevel=2) + return {'R': R, 'specific': spec, 'time': time, + 'bark': np.arange(1, _CHNO + 1) / 2, + 'Rmean': float(np.mean(R[_skip_index(time, time_skip):]))} diff --git a/tests/golden/loudness_sqat.mat b/tests/golden/loudness_sqat.mat new file mode 100644 index 0000000000000000000000000000000000000000..3fc7e6584eaf92f4e77583ff5847c19930b165ca GIT binary patch literal 40814 zcma%iV{|4y(C^l^ZQE|Qw(WLn+wRuZ*4Wy{Q+I3Iw(-<@iu-@>d+zr;pEAkJ`Asq> zb25`eQB*@gRE(61nUz#gRD;RV#@>RFRK?!Z-O|a`f&afgaZzqgMpAKC3sZLsb5bWq zeo}Q)cT#0DcTx^EQg$|e4sL!&4ARrsa@74=BWf=Jyd0TO19>83Fpg%{j$TSKie|DTY>GeKh)gaq}G*#$N&g!H-TRqxIhq z*i_o8O=(_ao1ySB7@2;VI%0ZXMK3=D5~dUIzV8H=G+(-H8m zSCZBT|NG+e+0je)hA5+X%ek5mh`R%PAaUv(Y^3&oZj#*w>rE-`|Wx>?7}kYSOx7p zO$_QouB_!)J7!_cw=V?xZ1fI@kYudOGx6;=lG;RnYBcy^Zv;VfS%rN1Hu#-v;)HB- ztnNTw<*z^Q#t1INZ&+>sTLM~7B#+%e_zC&S53c1Kn}#c`u!_0%v%ElmCa*IIP+M%` zMmEj!`a?-ITtLNHQ^75bo1=nQPS09ND?ZJ3&K1)KNvYN7@|u8_#Y)Fn`kL~CZ-vb# z#de_NV((e)e&91iUY&G5EpFTo4JWI(X9ksz$@n8OIt*;smwS`t(5Mo?({!mhPXacUeQ(eX}?*_O`xkUSk((k zJZZTRusJbAe=NqyoKCOuW4nFKmvT}ndqX)>QvSNW`p0pMwrxF%M}idiBv>jC;B2Y( z%yn+%UhHM4N5D%26mPE0UQnoH=3XE8V@Ee@k z(FIoqZuBh+)U=d5AupdXFkj@Li#x3o@HDXX1i1|cRs&OQMU570;_+1rcjfEW0Z!j3 zYo>0?U+d+Q1d{azZQqz(&zHM+8CMyQb;6n}?(@P^9>*Uq>o`|XF&{=oIvz9m1tQ~8 z2QRw3pZ2LOf!AG&&vn4SWtY~xCgUQX%Psfy)MtfCfnyB-F$_QeT?#}%`*ke&KU+O< zX*bG?Y$134YvlIEX-=Pn(IDcNh}rnNS*X)m>(~m=ipDYhOeRodQ0C^s)~HqeIc?$F ze3+P1IOr$a^;h5SxynDYs(G3l+!m-H9~pKb#oxtWVIveV=bg7bV8E#gX@HY(uD1u+ z^;uKmtOJC$a1+db!SDO~Dg%!5k3PYu@`$)tS#AiW#5O~b!udKK{e0=T0(&YuA^CPL zG>_~AmPNWsC8ijAIb}1U(wJbxRC1&8#I@3H&+cd~Fwkz>sh}#hhldjmvMV zr6%%6w~D|3cs4bC<_hjWt`)jTqXbi@w~M${Ou(a})2A)uuwt&`lW{_Q+pF=0#q#s; z+xDGaX>e)n-+)W8Vrg5BINPhCf?Loft?LZntWpev)F!0eWoE=@@yIjR z1fvvM>u7EbYeNIVr-@EJ=b370@~EK|Ml+lCU0hczq7@2^Gw`Pq?C4}Qdla-XXy zUVD4H=ey|mWz)sP2|^AW+k)jzKFO3rEmbL!TywUuS)8zISaENj7w{EGm{;&VN4U1VBT8aby{ zpT(4t!IYAT>CHv)K>2>|lN@E+KEW7-8#wy%@KBbedYHMp0#f$Je|R+gLh*O!n5p$T z(RR8a3ab5NYRSI}+BbV|5hgGf9e#W7;b~hiVTLMxdrzIX7;Yfio!P8cTmg^G6>9zh z#RXjfzyA~nKKNM~Elxme7Mm8jNAI1e990^#0w>D6OBvk$_Qi;u)}+8X$|q|Z)9q_u z2;i>)_BBqqJuy$owTY`{BNFI1!gOOBiyR2CR@f;RqpHaaBAC6{VvY}8h44n>J`ynv zO<{kDy_Y3-zsl>FkIr*>^2MNiCNb$wrcVTiqI8w{VJ#MehH_QjE&0>5gFR>x=ZdtM z@nNVlZB%=*S*^Npk+PYsc$znQVc(*-*GCPq-|UFE+HLf|h&RMWkl%81BfzPG#{E*% zQCI2M)HXQ50)-_L-vdVbB{voVDi?6iR_7cId3fP-#`o8Rk}3X%_n{eEgRuP#;%F+Q zBlcynL_3xPVMj`<{ojdTL;=C0>3Z_Xzf)scQ=QnTBZlU>AnM@0E6 zIalCqzxW0zA6EZ~>-y6oPWhZECxca@CZez;@d zVvYlMl`Fywr2G}+#h?TU{Tjq4KcNVHQvVrS>LOQX-=H|`8bY4~5dD=&m_n9&`QeqS zByVXm``jEf!C+O&&Dre>22lt-tTfLb3Z;ee<2Qv4wy997PEgZptc9p?!=1WKvGAJ zSr`H9eyB<%S-B;Os9+1t6FA&L0Lh!}ycqoZ)pKi)vSuZ5XSZu9uKkH?JizrR5NE^` zot~1ZCd&8ewM@-rK#=V+LY8z!$B1y}B2^VVjUh^8r3tc*O#ePGj5&q<0fvmHtJ>ZF zK0ZM)C<%$StM47HhHWFv1xt~2gJTm*$EoC>2L`AoiJ*uMw^!7wf)kiLLPFUTZYg&3 zf`_p$S>ku&U8$adFFJpGmh{g)d1!CNO9lo$W=(}pC6x5{dWaTZjaLn2 zg2XtKcMW=cfXxM&Y{o%pH?{F#M zYc`}M;@PW86KZ!IgWl3>HUW6uX?(|yaE$FC*D^#Om%>i%`6lRvEt&G^>F+`k9{hFm zyg9w;_#yAZZdwDuipN5_uQg<1N(&rXn!SnN8_jMWerDb(yZ2JQ{q)2o1e8^+=`MKa zYlo+?$@&q3Qj94__UUAJ8Baql3yN% zRpjWA&&GJpIrJlBbz;wxUqI)Izceb>MWn>zdYK}0>(SD{mk6n~l-r7c>4)UsWZn{u z&3}=`hl}0eh7C@b1)ItO8aX23gS8MO`E+D4Pyziq!DZZFvfaeC+|CrFCv}~!1uEb% zX*@NylazmJe>=HYB0)-&C$;5V3xUyjoEii7D2j@(y31iOqh#{=c=g00qq1p!gPi!$ zREMYU#3Ono+~2Z?&L92!D^6+x(t>AP({!U8ea@Msc^>K$1UXNn`T0-UtC-0riX9M6 zRY*y(gzklGVUdx-pzPN=b8(`*Ep86SPefnj)g-#MI_7PcW#V$la29y!u0=|Qqxmn? z*52%S=N7U1FKh}*qts+3j&glwuN6JAez#o7jE8CeTV+#_PY@rA8kQ5mveN+zHUGc^ zU1F+oKn#I!?tyzGI)b%{-yKO?NEQiuXye`Zy6NP(t_a2;lW^p#B$mNMry<}8zD4Yz zy{0;{Z<9Fz&pyaoLh`r5Y#ul1-CD9PWfyEQid?zBdqrd)}LOlU}4XfAHXP1dZx|SH!F;~M`w{T{>9})pVK3j99 z=fD$b^yW@9(OwIAA9Cb%F!c5qr_-S)s{fy?j6ciJ<{4{pB>mj~$`_c6^6kR4ZFz#+ z12@5Y==mvZEtwivc%R$+JvbPNji>=)&3HpqZ;*O1**^`9zrv&SFwzbe%x2QL#hOkN zTq{50=UlJ~P*b*l07`dZ3F~VkFEZ);kIji*1C5)OB`c6;%#XMKtbf8+9OZRf6Sw}g z0`_R~R6~0f5Nw$686v?T4A}~#RgoR*X_$Av$a}wC^g4SpVU2Y{v_IHYIt#mikQHa? zP5#B?GT2myrT02v-~%+p4m~@nPh2ojH!vlCJ7dy1-^8iC)C50h?y@u&c!RBbl}&UL z8T~ju$U3&#N&8*D=@kkeIYE}b$;0y(fXH2V2sbKgL7JtQip`-4FRmzEZPOtyIu&}h zL2C%EWBc|=dXsIfxW ze$x^jeUX7QaD+(X@r#>xyCbPnk^*6i`SOKOC>Ji*|0!Fi4-u8*E&o_oia!%NU(ledq9Lx})SJrq(k%WLs6K$j%{^xqLqRJ2@gvirAe(M%hqLw9!!>8I2|D3HQ4Fu5P;i2rkVW?!1S*HyYE?>W9%a6b zQC9X3gf1A9d65slJLCV&7_@9uG__6=t&}6t}qnQFYw^HEq3C$r^ zcn%HXT|`=ps#gT2#As<2I(NL3s7u;nt~@1tG&+HWUr{<~$`Zi`tkDI~jbOG}jsgF0 z7OTi#i*2Fl!HGW*~ftB*;n6MD&~i`i!AfH>rbSy#%BLx-czoA_rbM_qJ3zv}q& zxfQdspiwAo80WB7RhW8_vDdgFspKhmppLe3If4jNOVFUjQx7M+ha)uz$3xv)Iq&dd zohIMuf?a9$M>v3?Z`)Kr!L+!sAhG}wMl+m`z3`B)0-eL^1Wy;)miqP0M`3Tz%4He- zVuH|p9rYqlAo_QnH)-S(70e(@Zo;(AM1KaAzM-_I>ztlN=~Ad6VeU@Q`-FMSZvpXc z?;U(XLQ$_g6@2Sih|EK{j&6{L&EllD4^4v;|vN@*xcz&K$G(H3? ziyIHOp)O;+^QSnsp$f`Rs#SFDJ(Cz<7{#x=03O3<`acFCF%C4nzRqTR)YDs~YR#s6 zmOZ?$2KqjHw#batmorv}Sl8$LA|&Q!J6|Lz=Z>Zs&=qvIfybc@{D{Ljs;kAsB)45e zZTocWFm!=^ifJYa$Nn&w)qkJJLec?^fZX?35!6lIHvgJqzCz^OgDY%{ba$uWDyq~?-0@ZHi$XSP4y{@mcRRJBLVFa#HirK(OBX|+!l)^$9TpVQ~Ff;zu5gVPuW=rc6COySEvXl@kKVK_VQ4F9Oy-7if zczt#a<-s-8{$;+iuPaq?-QoH2DuU<|b6j}vwu-2_N-YW;_Y{|KJpa{96~UXj+@Vsw z`$L}J!>qXpvElURs5ZC=1S5g&{hO>Eggakpy-JyNaLgghtntz1p)OS0v~!8@Oj2D; zP@68vlRSMM!1~-wg%dPi2k^xu^?$oJOeRL|>dO{7t|3-Z$z`K3%Oj@t+~y_AE|0W= z@fsrnu1ek<4enB(ULM@?ba5Gkx;JvThP)?Udfa(RM&{O)G4p#=O1C4Q{{@;@S7?+X z-(ce&Pa&_0{&8!*>=1BwRtnIu&x=OR-@0@4?nyx=PVdTX-w)w#8&iD1*;*P2t#z=w zSGGQk1`(7luh^QAYWt`I%D5A56>@RpZ4c!fOzj^^*E|)O(pEOK3V4+e{!k7{Yl(@S4F2|6kk0lGV)&9tK?yh-D%hB%l$v!O znCv_z*xny>W9fd{&VLYZ;CTzgM6N7Z6a3Ks@%D-NDwtQzl=H~(_RW~tB9ITa8IZ6P z{?D)rPAxVNk?QHx)W$c9=H*_}l?D7}O7Jeu?G(L(1PV>^q7nppOurC5NmnSsPBL+HAKDcXJ`EV6wm4lT7mw)lSa+d+IRTs;+QG#De< zNBWQWbR#V>Z|3lvGseUUNN z6ngW=Wh{M>8P`?u@tVHkV@-sAc5kn3ksXqMKdi3!V=b6vHxHM@Dd^unqF&|4UKA2F zBYF4w@K92*@MP>71w^1zKCXI*!H@dJ-xFq1_lCGLh47wd zZ{bE2;SR=Ed{7X+qau;BdEa^)D)3(uuP2zM5k4Nwh`aSa&<5Ko-tWam-Gw)qw0wt+ z3W6fV*1ggJyKDS$R<@%IH8d;cYz7}bv?`WlibzHIg$>>lE*F#ZY5-e}_mk|-kD|NW zV$qSZ*C~uQ`{Xlz=sUh4)MiV@bq==voXdzSMKEd)eS1Yk!j%l%&hjVPcqDA+Iqw!$ za&8o9S{o=S>ZlY>XNp{rl1K42v#r_*V_X@7^?A~coX-|~gjsExQhEB{nBBaB>oQ^- zAx*YyZ+1vweM#0)4c_ZxXA={p6O(h8s6APST?ijy2oq%xlZv5a0d~qDHZQIMLF=ER zMl4HK2|$WLn;4q|$XH)iSl<-KY*8y^Q7dx_c0uZX1G=LHNehs?)fNK7q_0CT?nL1K zwe$6WYWEv18Tm~Ryap-gi$KAU~w4*c#pdR9H&hJ*Y- zry?KjbKFae9h&ZMsqU=mgC5dv>aF}K^+Y!hsM+R-Ea0W1nD{^uWJl+ScDB!*<%uvjKnd`D3iwl{nLkQur+ev5jeaB2oooR7fiH?1KafB4~J-mbhrd@1uTU((4PP8gjprra4Q(A*4`Y&LSa zcF4b4ChU({FZ7jXd&j6Dn1Z1y3~@)MrtT6Qg(!sPfcn<>u&OsD&nNPpUajKj2P zX}E<9J1EyZ+2;sp!>iL(ThLmK)hO2h4gM zx(VN)C8z%Ymw2B~9Zkn*^Y5$lkUh)Pr7$u?cmASlG-k^)R$NtG3C zNMxK-U%o-`Tg|>aOcGn=8}Ib>1o{N}%oI;CEj86Qf%PQIBKVm@S1Oc7JJR{+Y7E*Z zCP5S#(?zRNr@kHNKSZDY^GtCD7_@O9c%>iCpt|6JD9ccl1 zh-m>b?HX7{GJG09A55J~mj9~@*CTF5$#LQ?r-#HZ+!On0r6Llk=KBL*ZfjQ7b!V=7 ziBeivg(#sRRSch7Q+0sf)K*o!WDUdb7m(&LE{0&SoL)xnn7(YQ|6O;45Ienbuju+s z&jyc1MW}2uDTNdXql}S&K5eyuJNJXeOpP<&9yaqicJ;s6J)MV%sn@mx?Ghaw#l$_g+AZj zhC+=S5pJ4y#IQQ!+9$i|_w{S{XO1IF)Px17@rq?NAkr|!#&7V!qgK*a=x_alJ(QLw z0>|OuF+q#a;o-odT5M-O2mc3V`qBfMQN5!Q4*>(f)&@5SEZ-xMy3Hvh zZf0xEuAxq!%D3GN@2nzsvXcI+gu6zUYv*A}G3NWe=)LJ~rr4KjT-}p}^*YoA(A3kI z*qPXq2M zgina+Lf?W`D%?`nfjezxJ@8C!d--s%x^jP0Ce@?#!Ust#(#G_V<$noNZyC@II)>lf zg!8XeC_7Df=@%4VbLFrbJ$AJQUyKKouJHmrYybEiTKo0%#cw~UR^gqvSY8j&C*qo0 zI`PYHrV39c*Ui5EupHnqYXB2E2R@+}+4xO_3j^fv6iqnJz1h>379G_A2Z%ioIJR`* zR8O9GlrBW}1clvg*}tit5(<(~&kfOk<`nwyq%YZdY`)A56@zU6AsZV@@&`+w4)Y&Z z21ig3Z5S&@4=%gxqiTOgW{=G1s5!k$gmZ>?GFi7P^~U#P?@`Gh3BYe`=Iv53uSV_! zzk>bo1Jd@6Fm}D%-!0)4mo<6@)ewL7>TS+ns;eMLGY6iJ=sCF5Ed=68ZDJ5t>oEQ5 z4oMWu=i=rwFD`Xbqjjh^50Kg6Q5QDPaCd1a*#R>AJ`A<$g=SGY{m&LFbL^JG2T045pW^E^Gc^4o!vJ-%`P!kQrd&7xH0e5aQ&B?_45oSlAK^=jSP~33p((f0 zQL7IO&%fa9Lq{EA&X1vbk|WYAv$hWe2QsMBn92&}QPq|e#V7V>6WND8G1Rc$-)^FC zhoVNOOg5W-DM*oe=_0kJdQkluxyMK=P7I$SE*q}ck&TzufAt)87MShEOO>o+lOz@Rwf3Kb&>3_mX^5DEwJrSJ7%KrxE{n-(-~PM=B$)rr`!U0a^`R~Z zvJpp4v7yAMdh&7o&5F-tvT1c7#djR0j$y1m02cjY)m2Ar_Gus&%xfGbY4@1?zrl$b zcT*)|giuTo(3X0-DrpSDzEHF@i&F-;Z9~{+29|J9?ky@6Xi=@N6mwWukS7Y&w0N_)5aRU z^ijze2qHe~cl!RfI+l#Yn&v||N@B;y$wXJ(GM=}?T*l2t_8C-nI5(awPs^RU>{!!n zl>3lQT^&ck5?QEbsLeBjWjK~hC*`4SBYsFVF5gmO*ioZ;^K7b$vzT`}l9_YL^N<64 z`Sv>jy@@=4%+`*=x!v)R)=q;TO6w6_uA|-kG&`MHKHAA2-<=%wOnUyYZdl9{X#^Y| zW4vcuA2s;XcCp$CWrIiNnPBoEwB?2!AXXE#qaTl}F2ficCj8!KZ6?t!>W&eBvV0h5 zHrb2#2l*d6h_{KN^=2XF;;KSfC*d>Z_LqK(g9v=<{Lg=?JqD?4t8;uvT0sUSp$w`! zv%kv{24>a&K842${U_|lkiDQhz&6H6AUFJnR{3-@%Ch_^dJi{aJL*5-aBw#}rY`Tl z1+i;!4Ex{BeD2>3Xe6&Wb=xC<3TvCyY&IROqhx1dXH;Ppe_`kqN>aIwt@yjn(L?H% zvD}?vT_}#r+t;}U_47iahlC>hho>_~NEV9^MVO!{@E=!uF&ci3g8cJqu1NYnTki(j zR5l!52CHgwx;9JX(oThRC?5e^NXzMb3ZHs_@3;(GR`b%F8l`)gow$2J$SSOfJ0+9r zuTSbvH4E5nt$kmCO*mw@7G=mC^@w@&HI8HJ8l05l`63&NQ=|PYkPoIvXWjlhdhU2Y{u2;@hmr{ za%=|N2IPCYR&3ne(VJg9W0&hFu#Y%mfaeu2^B@Gh124OEBNWovQZ3l{&3Yrn17gLAOh1Imd1Y^=y$a=w>*T+-K|2}6 zGTZd_zI-5?lLox0eZOe>tepTeCqZ?-7)GlP7T1X;A4_N>9@EY@-+2F;OlEE2ZGVf! zRjBNifD`XtJ&z0jzOAme)T%Xau3zg#tMsbQJI7>9caD6ojKXXcZ+=Lu++X`ug{$C1 zlP^R_sVyPjPg@j)DadnfuVhc5o;0i%e_TcEuON59OD(wynD>@^Mefl(M709`Ex+!S zqgojF&6(WN9kC?&YlGrZBfi7j|Ctp4J$nlV)S5+n|1v9G%Qx;v7GFnZR0uC!2UDoEW2EC-j zd~gZ-Y~Rc0=Lbw*^Npe3^}}0k&r8Q;4QS_+_j2c9!Zl2$OT&U0gHMcidh|8kn3&Am zIQI(ni&g9~mxDF>l!(f}Y3rKDh(qp72PTpb4v*?8^Wf&bM7$y8uY(CTWY8?4VV} z{De~}45sI%orFUv^NSw+4=o{07{FiH3KW=|Y2I8}(x8kgc(Dbeji-!yNRlzh_^O1u z-jVLIny7^O;_yDByOTqT{&nV1^TsCXtRH7lC98m5*hMaNSVIH1>yd7kAEAhz=V&iz za7rCrT}>IW%BhHMWPWYfkC7VI*4AQPw}IIh%{2Q#77zowpWN_;B_g?%AVR{T&n_v0 zVpChxqKpO#WHC@U&ZfRdbRu?O%c4#JdaI{>=x(uh+60YovVHehzc?es5@b6?ofS|| zY;_HElv+k1RG!8#*jvHIVwHd0B;=9_dKoM^vt*Lx-&r#(>RCpCc-9;p^76<+k$8?y zFQ?E~*jH9Zgsk9bNfcvcvZ%ZtK~iB1R?!khu9o_{jAXkfQ(rpv3A_wCoO7l^(u8T) zuZ$#o(wxzrk;EOW(jV(OrKL}n2__;w`>(V%`4!*&#SgRA`Lyl69+5HMg2|(mU)*o6 zf@pFteuYb{%HG@S4@W=|>~uuC*$QUV@Pi2g##JWN1wX@^ZDk|U&`l+S4a1(F_B=UA z2qZeBUpDoeUKo_7S;ZS63X15wOE~lq{|-*mN5j06Z-2fS3@y;2qxS{01g0v9oj_$h zUWlRplKe&%@s$w^@)iJ|(V!=+7)7}d(iF~~FN}^3$j?n2*&;3(hLViNh~RFph!U0V zWv6_jhcVv;cJVvH@hTU*VdYF_e=qnQNxz}D;2JmVt~NpX6=Y}B6X}3IPfF&rcC4|s z_Z{+RtlfUp?j%ymb8|j*dL;TZ?rUCkIjncBJ6kmAZ;0d|rT-qliFaxg)l#fkexxF} zaY$one@I2zKNWZSHnoD;bvTqBh`CD2tDRWku#%2)j}QqU5_pW6QJ{%z9xS93Eq#I& z@FX-i)72eIn6*NzapGw>MEw02b)Q8f+v&&NUbfz!Oky}k*B`=|zkTf#$A{W(j|W82 zM%M7!CH%8H6B1SdG*TTW^e%k$usK&Kb{69YcKI;HV_>5dCTc~x7QM4slwHEoJg;z-VHqEWbuS2?hnDAhbKFkd z_46p55&7-@8|)8X@bb?FX|$96tV+>pH~=@>rfB|(H=u3LtbOCkN3be-p}9M9EBx&I za02QvGtREft}Kp2%_~~;>NXpoW}otrP|e2YhOEWYVQ*C!d6+0R^YE}Q8TAbbI}WXF zcfCCj@<=^Jze;;-we$Y{cvs}%p{`XzJ&^ZyUJ2Yxr=S-c3sm1QwVi~hGe_oF|4A&k z<)b!m#ERlNV{|H_Gxd86R6hA;bpMx^{p!~uD@qlepQkYT@$j?DTW8~+?uMRbea}T- zqmQTZ1)rBY@}{r4Z5c2YLSd&pHJc;JWK(~2EfE&SutwcY`1?Y;oim^uhVa;Xpf_MS z6=hia^hx0P&a`e#fb|c{SZVn%a=J)#P`7#u>nL+-B|+l;KGo=9;sVJ?M(EiB}*ynZxH1`XKu=ORe<} z^=Cb8=yr%KwVNGuW2!ji0RVM*v~Dl&_p(mXJ)xyJNt`gh`pk0eZc@Qfe1NSHPFPlo zNyxUeXk6z`mrJBNY}Q-)5iSN2yLuZwl*{g}Ht49;vusCSrN61OjF1#j?efHsIU z3s9S~EKSeh3sS>XjSIo`$BEU_jpLV^Z1`^JQ=mmxWKF=OO&&jW*q2o(bCKOWy#G@mE z?_vh^fzAnoW_vDm0_`W{#3HH->&in0u`HH!#tL0^SoOgiy1jhYy$qoxP}0ZJZ-`At z$0D7RC1F7*-HPd$YAzX2Yp`kn?#D-NnJYx~*HV?oB`^arW~OgtiJvZ{|91!YEcV=%>B&Iv@MN7i4 z%)YVHnMfUnnzQXqnw~a8__@`n(WAS(>EdlR$FThZEVB$+l}bj$o55+G4U7$LZgi47 zV)Doa{WbQe0ycB^gXLB_^ogY4OkkSS@&q4FEuRgh_H#GBc6xt3#3Z@T7<o^+_kXR86%ObVV-{apDWn{U>IJP` zUVHdJTBu#mbpII*;K|6I!t;4OiDj%;)BHJ>Sn5yC3?Ic>v+ zRI-)F@5wo5v8XhbN?5S{cZ620m8$}~ZP$~k;bMSBG2;$;qcO{|(LS}|GO16;Q{yVy%$G~oV@>fL%g-rF zB(84d*4S6G;s~4_<+4Oj^?u?n|NXX4TiycD#qN!pE$YYnW^F|WA47( z7eRk{zTs~DD!T7tjcvHjOGyVJC}yYDmV4=6u;(KJRu%M=mA9}l(sFQ0zf#{59zHLq z=lLd^9qH1;KJs}LD8JZXrgYr46eEuqckNodv2|cz?7^y9J#)wQ4oQN((uvww+u?d{A z18pL4k2kdvfjC!(U?5J*yb1bd^IJ4Gu?EpR9sbftZG9xkMeJFylc)CW8a+|V4Ez(b z_vKIHQ`PcIQKFjO<`BSXHu6FfQfal(P7J_q7`crS!Mj{^bSd2XsViG8X#7#$l6W@x@e@`4_KT+_x zHV7N>d|?IDr0ptJ1}zFYe-rok-ijQYwn?dwJQTvU^Dh|y1>LR<(j^JP%bU^v_6mjeSu*vn*Kt?=!iVV^u2#xgyd<-dzPl@ovjA62h z@!C9LD|t}73YEHJ*K+;p2}(@7^5hY_tTnmpRiU7avhKE*qar?59|H^8aq z6g%zvU?^pYzddL+PTJ(5u&-;A!OG_tKP1N${31g+ z|Gw?(?--NLNRp|!3<7>bvs^Nx>QPb>D>x}NE|<{nrFl2cCW>L0n^W~ZMz$L!HrxSA^`%wjx3QDpB_OB zUb=&TTi3}I1G})%*79>reenSIdz)`*yR8%7FIxWze1;n(9*301ejd6E2AXO6Boahb zIurtLcdE5V>GhvidcnX2Mgsc@!F$HaDd?3+zw%YCb(43RamM2uB!8bWcP!bqvw8#9 z&(s>GgQphd=h_BIiwoT0>|Y390|dj(A#JqNU$*psU2LZUOj_gd#_*?(^*=ar4tLA5 zg~kVsOj(YqvaeW^Wv=JJ^CV=Tgy42c!W%@FfF{`{l+(Xmm39?>0EMJ9_ZPk31 zk>xn(G_y)A&kAR+n!~7C4|4ESETd?uSp~b&XB4K;sB|y zgUgq1&n)#cLNNrHb=9Vw8*M*WTX`4(<=SldbCCf`ji}-4jbic99wrJ)nZ^x=_h){K zJ@a!Xgn3{P32}w{zscNWeMcl@*zv2GPS66F(uV3b`!(_!kIE-i=PXzw6``>avcMhy zUsuv%1WeX4PlUXK@J^vyxX5S-jiYzMfU~H zfu^RH)Q zljJ3%6imOPU4c7Q^R4R^ItsIIjJFgT&29Lsur&X+q1no9qQmb!3eyd6XiZk5arf^j z=hHh5V>+IV{1$Tm@5nxpw{J;>kSMGIRz&V*jBoY?ID|62ZllVG_S0n(X2&& z3qIQ?9FO)gH>GdlH4BbsI7Ra|&zy2T#pH6(djJcBN^jzl81MiIxcr0ui^Vao%QLL6 zS6`e1BY0Xbyy0!_{KbEnn}!w#@t8LTo|wO8aP836EU~T-Q!a)_H0Nh{4Ow56o%%Q^ ziJh4=?IwRMpGoQFpaznCRrU6n&yLnymK^F4z8qLx<>XvLHphUqns_z6hdW5( zBx6P%bmM90ky9!lwP=4N?01A@*%XPcE$eOV0V+Gbq`>$q{QKCRf7 z=MF=~@{dU*VGL2Pd;9|@GZDg0%*k(9+W$U^UDL}7Yqv0udv4O*z}L>d)bB96NeXq8 zENN_FeGv719p|)%P7j{2*OxMA`OgP7ZT4f)5c&3i@am!t1yBzz4`W6hey2L<8}zaRbyT*S=U&|9HwyAB_> zH)5`X7t^*q;-*4cFlx1_4gA9}JQ|DI&4rAip zr1BV<;nUmPWHw_|SdF~e+~pF$FD!-&U6J6?J%$(tIOZ&p+opUjqy*qs9D;OMVb{?>q5UGX=NJ)}OEZp-0;aa(pLlvh0?O%vWP3}8A zre+A?YeIJs92x$SR_bgH0y4ovIMUot9N;fvdiwa6RBz}03isdRLHQhdv1_%cH3)&y z3LR;n0<`F4>38YI`nZhAL6Mvfs%#5FtdIiGChfR}!P8;SRe@~ylrg}n>%M|k_bxBXB`an*SARQp zqC)54arJa`SpzaIVdjxr#l_dq5$RN(-u|_+lb=M#NAN`O%(&upg0*l6SuMFaPrc1L zE~b)qrj-)^lQ?E(wqs^yjyX=u%*@Qp7&9|7Gcz+YGc%uK z`}o+K|K6&7*t@&cdQfX-s-`rbTHQS@sm-~ZGIG)%YW&1Z?^(ut>~`AHiL1~ca0A6b zp8>Ien6L78wMCnOh3-SRbWnEK_n_V6P`8J23M>RXT4~nH0e*9C1;JF;^~}(~IHQBy zcZ)Lqw}0qVHUEvp zIK3b6cv6=zf1k=sK6%IiZ{Vo4^6g~KJ{UZBQozRkoqXelkRgv`bC7iGjh1dTEyrNR zWu*tFQxCXPM3PR{x7?JVQai`MuQ~W&mHWP>B^Q)>u*woi#lGS}_?f&Z9b_7;aCl~g z_vXY2<78Sim2r~>$%@9CjaS4AuRO~9`Z#_}4ej%O_K`WP zBIwZR*KB0B_K_dD>Me2T>mM=8gcE&Kv;3Nb#1actZ&N(Q9o_4{M&TF-A87-?(G55%LBJ)-n zm#y%}p&E-CpZH4CEs5u4L`;V}2O^TYrgAhwynjkfM{^Kw=s$T328q=k4% zp-$8mTINW0^b@Z-vO~fnSY`~Iqe@c7xy9%AkDVz;9FPeG0SnlKhv#)ZudfA*!n{d7 z&q~@ziw;UJ+1MP}F`}Gl=0Ld=ElK`(N$3s#CQ5~~FW90=4;j?5xX|6!&;oj{I zI9Ky5kAM%0<2@xt$I82dXXBjO?F}q;!$z1ntHZ22|3}h5HHY_@ijVX*YRCiC^y!27mN(Sr9)`(P3I*}xl!~u*#yW@P8-Zp%+)z+b!S0YPk(V= zm4JGoeG%G)*VlLZ01+*?rra#QVudk79pvQF&qNp+n0P*&3LxNw?XawD&idEdD#=-+ zK_Cr%6s49UboG#qBngaUFnP5mXH55n@&!9QSwP*uo>G|}v+2()8@htgAc7E0#*rg& z>LTen=T#@ZY$l}L^*cgkpbJ$|s<6pdsL7?4nY^C2_$*rJ=Ox!F4P4Mb2rG;$5i#MQ za9s9a07AdUKwn&hZQ_h!_<0@;P(g{b#yl-QG1vY1qkkP{hp^SX>g{XnA5dnYXAv)i zz@2f^Yr;Kpa1UKjzK9?u))^Mc#QY;G!=55CC_Kdv$}+i2dQDR^))|I%|8-r(Hk=zh zoG^E)*(B|Slk86LiHQ6({2T}80#rAmNkKw1!l`^Fm>PBS@gKE5sOy{RVg9BniJE}ItF|4P}_ zq>UU)&SoHJ4#1a8kDJP`8sncCiUqa$^c*o8DB{%{4Qo&Ws&>)RJOjc356?3KYUydo zx!Nrq!N*(=dB7nC${yV$K2T`_Y}W3cSoiEuu1>O_Ex)m07v89NE~~-;Ye_E%5*-8B zVWwl32BVGNvHk&p9!hu1uytiTATa{y8c<8Iz0}Ao6iTBeVgn0+BAbyjw;>Ha)8c%+ zZA5p!PuBib@)#tzY@gow0HsVLel@qrziT5Y;J>qdROBjhxyf?Ia0>F2^4F+u!ngeO z5xR@HCBIsL-ln{axX!uCXl~I?qFYuK^o#^yF}jAFdTg6^@x#?Uy(8aznC0^)JXhOz zn3uSE{Zqp>ozeAn+$D|!D6SFT6YB(-ps>B))2(s&+{HG}T7cSnd=e3{ytU&9lw z(5m9w)<^E{Gx@?*b8mG`=KJ5zL3=}9Ul6zqoh6b#;4dZ1qGn{@{bPb9b$=5ce@bvb)hOWQ4_OxM*sPB*)^t;e4mq<(H{ zJQpIK9XV85{)wflxzs5=hU;4HK-kt$=9PR0r#^T=uf4=X`XKwsRm1Rxfqu<^KD9%) ziVgk^gNx)Q41QP;-(?+ediIOv!UEf&F_x68`tI5?p0tU-7ksXISh2|QIFKKWPBr%) zN4C?#Y-32Jhz}d~Y>GWXc22~%Z$|RmqV!iOFAyY0deFi4Ag2zjKBzF9)R0DeHe%N3 zPU*?K|0fU5zf=&*UR4Apgd_LcC7dBRM~Zc-;$p#QqLvij?*Vz#v>1@TGe+2(NebLW z#g-z}MPXP*aH^8A#D6Lup^4y${5pT$^1L~9wmZGLS>JYYb#*l}b9FTpu8P1s?6T{n z92!b{$sO+*FR@?+&Gg*fzVe@Da4~w##?E*s^u85qFjYA1za`Ih+4KtDy}!6233ccR zeVnzZp1ecQ&VHB_f?mGrbRB=>A6g+_cm;Ue?4AZ=+>{C^j`K*ZF&Z*-Nm@tV};lqclqLqzV#bj=D&zF->&3kD+ znShjMJT1qo4LEexZ1fu``+b0gp^=Gbaq*jDjmXo* zjRx)*)Vyx-IP#RTIQeXzjGbon7{q|q2zas=Gq3kf$1b~K(ykfPs9GQ27O%wb4)EZ|&6j6CWt@=L3)W{L+T zeCXG_#-P12owEUkDF~9D%$%YQDfp9~|9V^>n`kGwG=l?7ypu5TkwBB&G*+$Qpm0`Q ztS3~jILIyfF*b+odjDL=^%q3?J#%4B#z|WDim@1A1zZhn!(m%(U~RU?aE#aqC3eMD!)jtx^?L~n&@$#lg*U4`8(t8n#xeM#$PWX*m{F1YT zBWuOEg>Z~i7u#iQNPm&k7TYBfcQZ%)o$=y`gXzFs0id};zwaRhz|l&v?hxOQx z#fsW-AC2V;5j-5tks-A^o}q_%^E^d}&Jer}lh`2%EdFpqd@B+Dz#tGYitSj9Y{PFa4-x6gW=MO6_qRul-}!ezaEHtNn1bV8qz&dd`;W z-T9R9XGZ63n&MW*vqtJB^GPT1o$I`t5-H1-M+M)S!GUXH=DDWnq1XJbGuQrBj`O1#(8{qal#*^S8Oa(}kaO@P4j zsg_Es(2L)5ZLnC9*z4u#%mlfV2)s37m6>CE*dn9XdcXNk+uXK`yEgo-=l7g7n@=A2 zyEdNP(N4stf-^o@Aa0<;&EeK)d+h%ZJbpR-Is=gf&5dkNw8!7);r~AcNHcugw|iYd zkN?Hoxs?{?YMc^&d^p?h^m{;`5@xymhxz-6a!N8Inia;4>_BwDKlDEZ*rz&KZVo%c z9{(RE_z~t5b%rcUk~`L(XrF(;{p|g5-vK5 zic%L<{$NS1=bXOk>hPx1(Z?qqa-I)|Chd%rxGsBQx|tT)WW$I*<+MR$Fpe8%el^*^lK7pFS1x2N zL@=-So~Io71w4lg{`eqfC~IOLx&94}Wr=7_F8g$f+6aSCc># zIq z+BTo~?H_@GLe-_-RzI2in|6lmADU5XxzUy2GA3bbVshe@2}PM54-0EL6vH!uLt`Bv zrrWB>C51;{Nx^jS1|QCbVZHk}pJm9+8W9Jm%`_G0y~JHuUi zb)p&S=uk60W%|pE^;+9je62YI^u0d2ls7=P{8LTp@WD=KCK4 zhJBW(0^7V7>h#-dK07+np^Jl)gi_*6zv-v27Ni?wiPZy}_|_46yBQ z`m<&*x|bC2TPLsGS$$iHr96Kvp-;BeA^9d%LZ+Y z0_EELUVT3TM|o-=@n(UNSI9~?zekC)WXfKjj3}|gdcFpc|m&8q#$jd@J35?Ofg~xgenTeh7w`r8g`Y`f=cX2d5l+cBDcnkGc?4`_8S> z9Bx1LYAQD+%o}KOchQ*roC$k9%QVtcV@kf2c51yvF(NW&#%=fFVv4^nwCmUN>HYco z@d!Qz%>qA9`BTrr z4k+{eQ%Hstg)fKqm)g}2-xI)D6oo;hdCx~EdW#1LMNK?!t~0Y_(>ws#gJhGxUOX_Y ztO}?^xP(Z8tiNT_6cujZlnW{T;(>h{OhR+!h}r!ni~k1b&t%^APtHm+pPbCzHNUU6 ziXP-v|5k`T<3x*8sk9WXkVx^Abl_VU5mFauycx)c88`Q;)V<5qJaTU|`jQn|%Pp%* zFtSv(e#)8bfWSbzUII@0p2^<~45+ub(u|JgMH39Wg`WKr$HlPk9xaq;y-zFnY%dpU zQ%qTxN|dzR{zsI1Q%O5vP*#rFryitg4ZsS#$MyY)A=HBHjRv1#U21T7&)%M0s?aLQ zvp7%>oO+;dmGfl;@@gtf6;8Q^b)uKu@|*Q92Z_cOpc8$b2cfTOr`;qAo4=G5akWK@ z{50cX5_0}#*2T6FD=lrnbJSJb0-)Dj2FS`%(ffavf0J4?g?!hD%4JNKM=e_7$zIpy zM4g&SY|D|FN6jj;Se(-2LT_Tw8wzD+K|Rnu5n!9t;Pk#0ijlnJTo%3faTy5XMdo{+ zvvF?U28pqJj&uF)s$KNDU>3Iu`5i}_mE(H2rvP#SRPKAE4e~K~KV)uFQc9~F95vv2 zRF^ifL{suxRMs>$v%g%gTq5K;6&}^PwIJ(?^7(cTZkJoN50TzfN0XJYCxUPDdtqhR zKh0!;`oQ!Y7`DS0)>?`Ro*SwZ0f?O-?w-5W9Djtlym)#&Svwj-aF3K|l!|4_^OOLl zgBIW*#@S@CZ(T*Yy*JNPPKd^2CD-|Z*FoH!Vxw5K!uG^xQBtb}q6wiz&&BDmS3-{EdDt?6nW!Y z|Nd|H-M<~Nl>vD6TLSBmn*jO>N-yZUn%fM=>KPcQL+H=5>FB|^(;p_9zHkYUy|9`d zY-L@8Q<(tMuZ2Jm@d7*JVe}?q2XcaP@Bk!!i#|{$*fw@Rj|k5B`#^zyGtnhhAQj=G z)CNuo;_i(b@<}YZmHqq^gccDR>TpaM4NO)e=PHVB*`<`P-r8 z$W62wWS2t=XW3i{DPmL6iHp#uu%CNsCe`T#L95=ox)+CudD|r_q$`TBewW0 z*9nhIy!sDWWi@Z@@H5e| zS+asx;?`+*kq>uY#afLak!ND^-Zk%ZOow5q`Spz?=P`OHD}x!HW=OBOb`v7rNiS>} z{5Yp{g!g1-qH;`W1M_ZylS^_IFo32%E0~r~!D^b2Bey*))-OBi8PI>cNen8@Dn%MJ z5lOX(Dk6CkSJ*FZWuv)0CSBK^4 zaNv~zec7k?C)ILf0S0ng$N4Fq!Ci;y-PcczZ=eX@@w;QEycAfF!?Fk4;(&af>?)`Z zBSl@{mOFyEo-1F|1ifvAb1WI16%m@3VU%n4WH0x}Gm6luW{)pUVdf50< zFnvnKb_8S7cbOKWnwIiRF7KfFV?Ed+$2CE_3#4N*gTQHlGG8mizoOD7t$OkIp~Ydo zCuR_?#yW?GZQI`sI)Wz^`8MMWpJwzy*|JT)v*w-7yAx_3z;fOWkT?oaIq9g+5tw5lcz< z0cf_#6{Wj(AwxCosKb-`3wUNB*U2!YPR$Cr6#Fs_Rfga67(!={^}JjPh-R^k83I?R zxlP&@qpwxr;kIKKlp0<1mEEitj;*=ycE-}Ky4InW{dGWB+R!A5JVwhmLQn;@2J{rx zXk@_PoRh!Pdw>cL^U^`jUYzgU-7n6#!s4w_s>kIpbe&|ZMiHE-{Ov7&kYpBZUT{J+ zKPygCf3DK%MtKz@@FZUC5Us1oN5z>6^~nd($&%z(xfnN1UX7*P+4?)nT3$INQ;G@L zrrdb0g(kNO6SW*>pZ&QiE-||prj(mM*%yGgb7tA?Lg5b#ft(4s`3nn^(tMHA`X7R*O~2P^!U7adB@mZe{Za9@0dm?#D@sKaaTRz`GUn z$wn{K6uGmJm+j&rmLYN3Ow(&!2;zPmr-lcqB>^v}mxbfE^nRV?9ZCMovu}7&orgSD zU30psN_^M1L{&_EzeC6VO! znVRvBV8)c%8FffE9dx?_MMm(Xervl?+D*(=3uc(lAY7U#rrs!^bKdiTd>2BV>`kZ` z0rUi?d4e_5_!kgRj4Q&<#~vp)TR25BcN+lc=S;3dS{9^s7efY&_k>0B>nLsWiP6hM^$A^g&|0$avdc2HtZsT=r~ejg+=uWu%Rt03n5-1wB?H zhkoOv@Czs(WQ#r}E!Ii04^x=k7@rK{b)@1(^wh$IPG-Ni@jZol!QiJk^dXthI`4)y zKXFq6PAo$ay$!1T3v3hab)lfA2dLY=8_2Qvht-(-DiKslPao5bH!0^05*q;m{37FIc>p zemS0Lk==wJhwpx@M@+_wj(FQT8%92Vv=(p(xVH?THXN^hkybvF<{^qA?I-~G%AbCZ z(h^3xNR#LzqwCh5^{`*ZIlfe(Kg0sw@IkyP;|tOb7BZQQx4$ehPub5>AU{AU75!xh zvCpXH*Ar(P;hDMf-F{`IHf`8>T9X(kw6Z<*_-?s}*;k&VDms#Um&%k&RMI!Ba0EzC zH^qUp>S$+n)mM{$?AdPFcBfJmOY7Z0MoYf#dMPw;CW83|UKWU{9mV!A>?ql5%XG0T zoggS{OWoaAXAm}E+hLiT-jiL28otEE;Cy{d|K-0;qVL z4?r%)VLxiv64VI*SKPjKr&ZK%vEMpObYWTOgpzk8cw)Eow;z0~mk@FS+x;QfnxOG%|-omth=QpMeKx(2hwDTEHEZl)Jv?cmhH zbjH!WixE&!qQrt%cI^7EN?5q1o@Drq4)@7-GLJUzT^n=%5)=ZW)}4*NNk6_P^)W%1 z{1Ds~9AZ<4L{>2J2w~&c+lTlxE9508YX_L~#sxNeRK_+-F~n&F-fK%i4-g`EnZL|8 z3cBsgG_xS~<&qy9YYKeZ2!Xy2@j<^+wKJU|%#sLCe*#?fKDskX3^H9dHth2adR1pE zpCfZG5*Z#t7U`huN&tz{$9x;^dUZz;h>wJWx(zuLxms588|+ zP-{f5u6i)Mo?7kq(@Q#YZ+wa0AJj+(X4-2yJ2nINJ!8T@CpC&zx&_@ZdXy;u0rS&# zpR*k#)aKWE6}UK(!?jGBn1f1x4QeXvoHBz@$9uJuTs#IxxtUuV?<}B9F|sB-5*qHp zhNJ{3Q#yONd@JjuPV=*UbI(StLh63qlMY=o3fCmO?5w3}fY%h&3g%;sW zkrWbV!Iv<|KboY)2HVAeYlDX2%`i#ZRnUi*E)^&F_SKbxyZ6wCG;sK6mQv%sK> zr{OIEz-e_C7^A&*M}b*O`;j6O#89q0Im|_AtqurVFMx&#~}+ zo~hS)x=35`biC^{#-L0FkmK4KS2@E#`#@SmKQ{0U&WNE2Taboq$i z1M;aIPLohVG4E_>a-&5Qvz###3}ZDQqx+ys6M`O*fj!i^K|yubl2Gil@rjXJe*AQ1 z7+4BQ&WN=wP6L0u=NOq$G~vaFC3<+~cVczEa|EVC0uBMc~|u^Z(SGSAq&h@IiOUW8E|zO}0P zs4A4()630vkf+7qtPlHBq|>7`>NUFUjzNNv0jUC<9&>_8hBOq)Fxg{RUweEgm|Z*@ z7Ae~a=rJAP3;m=}Huin{KTbSlf8prkYt@XZ?z?vj&zq{NQb?{bGzY28jX#~X+S-3Q zE0@9Dj@m)d#S#{htvQOS*q5G>W0;jJS+be^Ta=01EY(9v?fu%+*0N`VB$eYxQ5#s$ zmpWcWXUc%h&s95acy|1`tV&gG`BNiD5cCOO`*1ie2XC>%ac zLVYr~gN>5u-z*WOrnS=hY34245G#NETTHy$YK8H}xXDZ5Z7EujM8)->&*sLN`G!09 z$0Hs5$ZI|xsL6>0!?OE2P=Y8J@P>9=i&fzBSt1h&OQGInt>i^38%N3op-Tzc|_PFg4PkCBdK@7vsI95>0ZT3Pv2P_dr~XnEFzR?AkO_z=*4d& z?F~s&yKq>S ztizvscOMD+8UB?=GgUN-XlXX2yR{c32ckuw@E^0zS^Hq>yvX^_9z|X56yWjY3rA*RrKM%J_uOFZG>rZ!$BSs z^bbp&Vi+Doo{r&UkYw3Ar&5I$`UG*Fo}WIHO9g@ir#u-DJ7(zY8Kz=&&EKWJ3($Sp z4O>cNKOg?plM6Fk&)UR9dHFPwkVk9aJ!I65WjQTo>Dmg#ZkwNhM8P@uKD{B{d5^TT zX{BdabQv;Kh!xMAD6CBl{AjE;_a2nADP;sY@WYvsOho@nzfNt8|9ViJHK)w1e$A-Q6OPRQr0U zR2|kBSWyp(5S{CnX7!{p5vAXe?~B(Ox(Auinc}nUB5PlZJmDJ51L_3T>m8u| zo~O63N|eJ6zSX9%`TQZTdY8nvkv3!z$f59i6Z;1E8C-0y#(IOD&G7p@B?5q}Js)UmD=s9(c&5kd_DqI(CayXS6xe=y&J{&00+Zsd8}%Y9EjR)XdI>-gaCCMPW|}#T8mv zqm*OC?3}!|3AoI^3|BPVUP=>;z0JdS{t8Dqul|+S+0bJ~8f|s%xM&|bngx6YN*r4STeE#n$w%_1L%tsIMv$Ofc^PCX+ zP!6)RPXPn9NcJ7bWfA5q3SZ{Z@Q8Tmrg$vD0*44kZF;S(=6( zwCFTxW0gjsoc9Y>rB;{B!1#ihpWV3-sH>mq-05%GySpUGZp86vVx+RLAq&*O?*{ zKZXwd{%d1;X!VtA8T@Z~B0?*QzMzi~qG(Yn42b+)#QsC%w-}sBRIqVO!}K7h4XdQ- zCcpRzi;`k~e_RoxjpD_MF7RQpdxjUtbIQB=42I9S1nFKj5-~lT^mvsox9{O`$o$z> z`9gS4J~$OFQdV!T21rja`0LTyKCwZ2=PLUE;t}lsz^Ao4X}`2D+c;bC|A{X90=ulr zXceXtdb7;eBrA^=D^~dzaBxA%nX+iabhi&OmX2CrX5@dR&7o)i7WxThU4Ha#hA^C8 zz?I}NBY=M0b7B3{N={nKRc=$LjfZ1>b;ad79y~%Kuz8EbW3KVWa|BW+YD=P!W1&>+ zU)Jpt=}-?8J(T_^AV|)0;iqpml#gCxzu)e`ZT^)mv|yjS;a@!UL%r7~Ksk%k0pJ3i zf@JftrnP~n)JTMFP~z*oi|00XyZB?OG~rSNe}@5}ac|{X_4IcnMx>Ie^O@ zeB>zB%~erEl<}JJu7HX`+}D^tgKrg`jr$ayk2RLwjG9r6NM)CpKm8c3*!#O+e%F^v zL}eY%6>N0zI`_a$-KMtR)2{eC^^U!zCqApc?i@~6ugbV%S2eFRkJT!TCu1Vr#?t54g!wPHEn)Y};id83d6XdK)Jue8yG8}F`j4g{QrjAlHK3c6B z@@QRid_A{R-(=Qt@%zMR5s~QBxW=pZdoBc+N!X}=B-%ud#mbI#!D1yP-as1AY%2FY zejzHWW82pLG_8+K!FSl1r~fuce_+*-Q23J?r}X}>eBZuaW&iJpu$h***J4SIgPjRRfPx z`vc1}&F&5#f3=v+AQVj#2`qulc@K{d#0s~(D8qI}^J{a_P+%ghPB!r02D(xl6+{LA zCeP%^pUW_k{BDIwssx1aZr9G;iIJc0O*L}a>I9lKiX96NXa^cs0!R9$ulQ zhBbJP6FEHYc{HDY2{MqQ?RXkWec*+yY<%L-l;J)S;ilNgl7IM?X<|+4K!n^)&ju%4 z)H(X;rRIyHU+vVp+P5|vg~xanuM{)>5BZKeg$|tZWJaKAhU>o$EZkh0)0JhBGIaF$M7p~DbO49jkWN3 z9axbGxn$ z=-w$dd2dq};fy|^dpr6!KT~zrgp}#KsNOB>W5-;~ys8vY$fJ06=d?xR*R{2MJ=t)5 zyiN>z(GBO+_FDke!UWDeGt#YMI}Gnd5uGAjCdJk~CColxNPe9oErXY&~SVh19NeiyS0=c!Vp zsp@(bw7qNLHSY^zN>wH-5RM+aW^OWK-_T4|H7-|>0j6tCx{D6?7=prWL#u>=UqPt! zum_xGaVoo3jr-pTGzX;0dE+uPqZ@9Wv7ZjIneBe!oZ2%)Z{=yJ(EZ)f9zHC}ZJ%uq z-phq(Ng*b%>PxBvzl$q6`mbP*gw|H-uB(F#ur9N$J1PHa$TZe??I3{@FJ8~nb>z{@ zgy!@fz3XR{wg(-zA%%>6g`BqUG=IQA3B@K?!-VK2W)+Q7vbs7E&ypAq?U~9Xbey$w zqx4++PF83Co!F$w5V5hSJD-~)YMC7Pn%Pm$+^GGOD}N9(V+?HJIJMmjdla1Y~B z^Rv#Z`zgtk18sQ5QM_rT;A0S^dN3h@!PU)xc4sbnEntxX@q6~?vJ{}*vB$1H)~D(t zEVBT)D(54<^CM8mv_;|9icruGoVSfYKNW8CS}M~K8ilhz&7I?U62c7#^+jbQM4ge2 zElgcKbA(gfW(A)(ivd&}qhgLI){DoNZLGvWIya44GgyK7ocZ^`D!Aq)sAie6Y5~*! z_8L=dP>NbP4_s8sqWo@i_RsBD>Ki-K=-Vr>J@T_3=R%%Xp!e+c?$Q>I0WtN5*tIaA z0(;>c!w1Ee!8bSfjx-U4v=b_9Z0S%;KYw-(dME_B$yH*~5t4A#32zQM)N1eJw=2Tf zxcF@a8W1adh46-JwM+Yt&-(Xk1|-@FDM>XJ+|-|!Qyd`jiLgdK{!s73*9<+fgzIZw{8u%a`R&2tWD~4*p)0b zs6a-ZIdxG>tBUDSfNU~8Zi;vvb1%)Hb1lL0e^rgs>ETzD2MvUx7Er`8i97Se^THmFWtpP^jpqKkH^I8(oZS7p0?h;IvQ zr=ffOuIA-V0|jVk>)3zoST^1k(^GbHvf{>p8o0-~W$`xPH|aa(_0wNA=5#K4a}ywE+VYTn8>s0KVb^L{(k7+} zT7s@e(`vDN={-|)7u1eq?dkrVvx?^c@~B5-DE6t3Mf=P#(Y)hM5e9eHH6Lq|v%V&6}sy;GPMBUqtj5(F_u&8vU?!6(?IC4}!bd2_K^`wPd4QK@bdmcbai|NYi1T${5Vn$h~p_nm&dT_w_s zIhR%K8lRHa&!^_3SoWdqO3IwWW{*~%F@$RkK3nK|CFG)xFs$0?Eyvz-3kjOBE1iR{hP_+;Q@wDFrF4tS_^)fv2;qr(#qC0w&pJDMe~l)`4+1D%6fmdf z&iC~DpKc2#ZgsJ!biB3O(7?Ldo~h8E7(Cf;N#?|Gq}#Lc+>$^x$Y2ioaSCl@Y<}CL zhR_7K*QN2*!&dAY>&hYDO7nf8ot@(g8mSN#BOo}#4gRKtIX&54RXB)Ywy~%ON)o!I zA-3E4CCrB$ao9o2LCa8J6l`QvZn}r;w4moT02we-Lo5F%q4x$Tf<^M3bLRKcQ$l$?tn>5&;=AEQVOcI>M!B2s9j*P4xr)9%V z#edzMJ;3Xzh949P3rxJp#r4>(zh<^Pe>T-N#M~)_IwMVN)O7boUH1<@H-YAT-B1LF zpJ-XPp&e}E2Jyf@mvE+E{e+)uJR8l_0_ zwZ$7Y&=fpgxDyMu+^qVy1!zO<(k2Zu$g1Z`HKUg=$YRsDFFq{vWD37T-mg>-E#6F0 zfdJ}1XkA_@bV;x=I+5{@qPzm`G4mUAw~2H33OyT;gZ(9dPT!x-d}9LAYWXH-XC2P; z?hk)0*DJRTnThpQjf?%UJ4Tg^)+T@ag-TvjZeAi?BsP_jM-g@ISfhw}^`7Gklq=J@ z+T_TuQ%0zA=TuDx-RZJ8)Va=@93HNQb7FK#`?7msticz9AO8=4FqI(&z=`q=d`Q{&2Br1ql3E^~cR zu^C3(>U_C^k9{Vrint={Fk!SzJV}my zGo|LE%}#^(uqA^JY1LRKY(%VOgjXI^wM{^Y||9kz41f{!Q7q%VB`&ivF3pRb?i|WSl_i zDoykVecffbVLOW=JqIR$W$VbLDO<2=Gu27?vWblx+q_+rc(3I zaTwq?k{1uEzlrV0BL)#|&%)pS8yy|ymg%$h=cPis zsGIrpexL^>P`z*o&WNwBZH@}bVSHHp`-D_XYO^!}tIv{wPi+iuy0dIBm2HDV{eM$* zuc5m)l6ke_#j)wikS@wGsASYS(d(_eCkieaxQ$lcB1ibjR zaT#;!U0p(uwQSm~Wo810sTBh9?bM!YWqrPRV<$$yZvBgh5;(U}FiGpmL*Vld8b(%< zP6v?d?)LH-TPL^+_D>ZjtBxPm1ZxlDWIgnsa<`3b)a3-^d9e65n0)4)SNRFt{WrLK zunnV#jeXT$?+|i@Sn@F8{5o>U_T_+11lGE)mqQwi+-4Nbzt0>f1!*_bWj9a2=RWn~ z#e*dS7tkC>5Ez>{^zz>NW-r)YYUB#ROwd=45cABO%y;h0N$0*XKZ*aL^=y~9?T-WZ zSKbPh{pNSvJnawHk22*OziK;XaDq6&Jt;dindBrLVLzGCyM=BaG}|WTTM{bF$jtqo z{{{2QBA5X2)#-3@Wa&WRe9`%kmIe{ec4lJ#w(---$0kfDi{Qe9r>*w9Ll~8IonH(5pFBFz(L14BTp&C$W8v1@y5ZW=jgBckWASb?=%N$>$=H2 z3s%+qg8h(2{6@Ati~3_8EZq;0BHUvgDW<-yg9sXp{k-%EvRyndd@Y?9bdXx319Kj6 zl~+8Yb&-&4+a_TX`WMLS?&vTWoLzrhpW;Xpz7t2rs|7`;pS`Dyoa#-KhyeV?JT4e*W_mwO(OcBw<5{Hu#T!+ zix6KxSP1gEiX%6){CH`k3blAFyrB82AGa>0&WjVelCkR>^R;t3;Fs!TaBo#?lirI~ zX>!a}gXg1%LNfjyjW2)aemfnE?cCvpqzH-VY4w5`$=$yjm$9@`Q+-@OgyjdcuEu)qH^O=J5-u=MLBNHw&$`RwW-Dch3ym{`CBK4(;} zT8QyE(4Cqq(p~C{IVqpGx!C|q)&1sdy`L7b$qm7m8Fl#?e4^o{bP%{-z7)jsV0Q_^ zq}gAj)^f@7`Ki9J){g~?u16eV^B9;14OH8?4$z;xcoWDj^Ak9U_IzSpW zWv#rwl8C8qL3>?sULFtWNxD}rBP0-hH z=5s3<0gj**A0a%f?(i5!+@_mN6r!mV%mfzmZJ@r@=6WFd$G^*k`WR7x|B8@H{Htya z`Tyzo^Zwi1{ZHuH+P^BsmpfnVja#nW!V08H_w#s?**!m!(Y$)a2G7gPnI+I@(-jkJ z(KDb4rHaPl^N-SCQPXL}&dH0FSwnC9ZgZS&x_8}rTyEWNHFy8DQyr_?T~H4w{bR62 zGdvvZ@xL1#Hnvda_0Uj}g2wx=IiTW`#`~?{?#`cx@x@cecW#??1MoFIv=>#Sof~Iw zE#u4PcGd39mg$FYDR?tG`qIu#wNuB}%gpZj#ipNd#nzQiyww)KSH+-1c$Gb@?e;Hy z2euDg2YR*bd+WAuFXqgPt6A0KCf^~U*4BA;W#cpFA1^vL&r<-h$JhGaf7PRr9|~UF zJoe|!TC-(WHjq-B%rs}x3uRW*U*%Tf-@{zU#jV~~Jr?|!Z*MP2!;`%Pj(Tr_vnNt} z+XX`0r@gFlSD*d-2D(p1J=2-M-ZU!7UBAwZu}_mMv`+)yk8N&)=ZYpjw;o%QT5p*b zZj>49o|oJcAh`cs5x0J~_7!!v)6^XUZ5fr$#v9Nae^%0s{HMcfY$YyuaXQCvuISEr zpQ4=PXGi1iFrL)F*Jnse$Xbeivv}R{P!C1@%?9fr5E6~M%YEnw6LyilUUg`vZzu-+ za*G?AF9~G27F8mFg~dz#9(VzrtvJ8BKHo2}m|dB1%MoCd$9{&WaRsLMgNH=P7dq;b z5^8BjY}3KVo{XZc(l|2fl={5vLGdu73i&zRvhbvB z-f(oB^bm}xxF=OKkM*yTkqKbYOo~?A-yfKGKI{ROdgM1{|HTGss=w-caJRlNL7I#| z=zmn`zKRv{ms!bZHH-cJ`b4QIc#~H2B>*T|WLFO-pRU)uyQA^-yAuASiyTc~$==W( zGIWd|a%P2%mj6mnF6S(;f)c;AMGJbkbHzfC=|mCz%Vi!ANP@586UB{k{7kw!6bg}< z07Tk$_lM}ly+%b5D2L1$s0!zcaWX4W(sTNo-(`TY0-?Td%K>Z~Y^ z`ZUs2H=b1mV`WqNw;D>k>k7VH1&>Z^j{Y?`)_ z0156c!9&pCy0}}A;I4~11b26LC&As_CD`JwAuR5&1poHUQ}rMGHFY0N^))@+SNFlp z+}-jhsv)5GLDY=_JY#1sxBfdC$?w%Rgl1{A?^!!Xj>yWsg{Vebl7O2&4@4Q|Q+Pl= z5BSP*f0d3y+1xxBQ{iF}oY!+)(oRM>=*pMuS=Br|N0{o|A^40*&;#0ZxAw;b)GpbkAny+3Rf`_qI5-A!d~!Ec z7bIwV9$Ls(YiCnfzH#*+C*VeJHdxdi4<%_QGWWjp);KeMBt4LkrsKz zF?m-YEZUC*sR?FO*S7A)wD4KH*56n!ZHGihj1DRd~ z&Hq9r^#uxLyr97Q!)`x!xW-&~trS^m6AF7Q8yT~}HzLXx)$DU5Nq3p`F$vt_-9!ZjkFXYX>Jn;M`F<&^N`)On-((i< zHW#U~EsYu$KwQ8>Dn+%k+!61rl(%N5dub@ZXowF;O ze|I2-*Va>;I3mBr+lFO#Svr9PNcE5C@5>ZR$YHz?=tXR{i7T9naVFVUI780G({Y;f zx`d@;`6@P-nfULa?1{FzAe*;3dWYMHO{RO3w8wZi7Agi{OkzqWtb*8*)Kvt2yK^L0 zto0FmKvVq4pzBOBb9sKm@LQ04!?8I!Ga}39`MVv0&U(J3@1j>tkBIrRA?s_@>x%U>5myS5Vj-QJ%dGg(-cpSi$CM$USvuT7Ii5cvj*lS6 zb6;(lSYARG4yv9W!=?5a$a+ z1-VYmYLr*xzVflusGVRYeR<45?@#8MkaN2{_1ZY4BU!HepjXQ`^n}pM?{2FP$Q?7f zR%Eqo-bARez%SnmOkuRmuaD+B!Y;ME!;S$a?;5@hqET#f1*uUodMNlhyvtW;Od{R? z%hs(XKWdFjVRE5Z_Tq+30lD`k?-5SNaG|T^vu#Mh$Yzc{L5qE;>zj7;SVsid>JaJB zo%e-D+h8d6ob@C1qav+NNY8?C7Sma;jAw(drtDk)I`WRkJ7f%m|I}*F288+mR|&22 zo>MTkTBG4qw=N4;|5&BzQ&6~*%O5MH;RJlQ)K_oLf&uICOieD~_NgzX{>nZ;E?`kJ z*UadFu|4LG^UD4J^EUL^uK~BOau_m-NQO{w{9OwC2!M_8!*;luGlcP4?Vi%^hjZL? zU{BS<5}C7C>Xb}@@#zHl<|KR!$Y9NW44!*eCkrYP{A>u*r7&8qWDWL$B;fgQM2@c-?QNlSYqmSrEg{_R=VPxHrMEW<=U;udI1J0A zmru|s4#W8 z6|_otpWe~L*u*u2ny{mUzdYY-{iv(knkhc9_zeC*HwV$?M{%faCYCqbZI5aEi4HQR zbdt}yA46gq8kNuv(cz#y|CG@JLs0_~=(tyU~^v`t?=Iz#eTxNbGp$mOi7@E?iIvFnlAm zSnH!bD%KG@7nF&{G-Qv_vvVC2S^@rsYw#jd(9k<2W8d}@I~@41FXp+CA89;lhp2cI9gqKe|Tpi z^h|$=^VQwz>#rdIK6``pf>QAXmpdzwL_Cjx&YDku&!Vqi$@}WA46QHK`Jt++9pHCa zs=JSyBnt=4wTZ7VD69J$n(4Eo%yGBUuEW9NmB}H6O5SB$0TCn)sZY>lGw;T6(r!?VBRLjk%lIC2pv~V!CHMC27)KFW{ z<`&bmt*+f8p3qsWR-~bK8?u}hrhNR-u-^Gos;&nCMZN;8* zZEM+RInPT#R;1_MvjOI>dadQ3U@p9?yBBEkpXzOfoaikQ z?nla{qw@K}RKNO3`ocUiSgC_PsGt45Nt*}eS(!{qx^7}mkAJ-tet8Rc&8C=GtC0Hj(S& z<9~(IPSC2O7rMpHB@jzmeA&*Tdy_aYq@`@%9&TR=opq~8>@e6Ja_rc0W4D$XlsJON za@~|Mv2Vtd-YPhePLwAIT~_P>F<+OIc)PhsaBi}O#U$vB=HmNjpOeThfe6r_GtVic z-*xzFH#{u!dc_HpOsVr;qYEhu@_|>b_raszy=vaE(;D&GCyId5WugYsw$BJc!$WR&i}PPx^Gw zUrqH!=NCgP=LF4?K=7&O)x_GI6xtGoDFIAb=}gX2ift;~Q_9>}%mnF*Hz|r1(KV!5 zwWPn(wIhhe3q2<&@d#fYAKE$n+s%T2Ml4gNEKbchh?fE28`xh*3j+}pGhJXT^VuLz z&R$Q6HXRdD@+lWUWbd3Gb$q_^O|CMH1j)gC7Xc=RiU~6C6=jTTwQDyv*u% zNwwsj3RCpR7?>o8&Z5~H zRB;`N&f?jxcf%a}+1KdJ|G?FxX=_P6{fG&OLrAO0fhh8MzDTw?@tx&6Q%#PK^f|9d z@BdkE*>B1pj|66@BAitsQ&|@;;x+Ytjwn-{)|dfQek|fHD=&D_5QU&jRB6aQ)g%F< zkkw_MYLf)<(~g{J&j#fhG_}IQwn6<=aGr0*%nlPKc?P zt;>klwz(t-L*;DI&Tfbsl$xLu_d4vKn6llO&z^L!85A1jna|#IOqkZ%(q7T6Wz*SW zrga{@zl-@{+gCn?yt3pa^uql+vP$q*NHHhvi&x8(E+MBvri~)AD2TwRi`%-d{J2KD zf*fli-R722#7Tn=%4V0r!Xl#DQGq8DFVZa_V4y}REeUYvcRfDy>$+HHFx!70`H9u! z#CPky-J+F?=go0E?*ek#c;T^qap2hls%Y7*KjCNedHjIP_}qb0v@Kv4+J^SheRHb8 zOKRA0H(iG3m)UFPd~?HH;LC)znu3OoRv<<4Jm%4UZC4<;n(P5CLmNk>csv8xz?W{l z{k1Md`=#EaDv(lz8PfNu3#D@JEWoy$SO9m^w^68V41I~wpKjZ}X|QOiqyHz1z@_5z zuzMTy+O>BTmbH9BKpLOdDBZ}l@0+XLPk6D!ngv*auey!!EZHi&ABR#J%T4-V*T%JR zw(|!1BwyV`qlq$dDF#O8z-)76p}bDcMXDpAG&9M9uVXIhjI;!=x*j!aa;fn?sE{qmF9f8pWs zY5QlGGQmN*u-Vo@pl7)NhE7&Sq63?L7wu9az4O$Y_spbFyobz^j7$PyI+m?{R7$9A zIPoXe*td_U0NupQc9H%$P4v21Clx%`izNaS0g~TGwHB{;8dqt;WuM!d+^u61RlEE+ zP7fP6E2;wLL?$?+{zhELUY7IZ)*Pndi_3#bSbKEKLt@XW$7_0#fI2RQ>q4D*c|-|R z7#sNtZ#&NH^)}r{8%<1`iCxfxV29I8Qcrw5ABJ=7YGalwtnId2qADJ9?{>IP+C(jN zG&Ua&nE@XPUs|pFCA}=auNyzwE?!`LfdAaVPq#Fu$suKKSY|$bpIKIh6L_Jd-mT2w zVtv1rjd@vHdxfv|t&KJR?@0MY@!Flx0nf2$S{bkTcS^WNvF;PZbp9ji*=a&WaoSvR z9mn78%$308#~q$;hyPAE9Cn~jcISK1unN8-HkIo-_S6CdrYB)WOjX3QUrl0?-HiBDs+f0PK9#|Tr;+fnuD zN`UgFoHyt9PkGH-_s_5gRgI^HfRhJ+jgZY*^Szg*vb@6LFG7jR@UO_+stTwG@co2r zL7nqV2>JAAw3`Ny+L`2uL9&(kBMJW*5(T)yc7Z8Zkf39`M~0~|`!*OHO13Ra3ezd- zN$}hIT!~|o@tI1N&htrLo%ih*=P5Z?5j^k%Ijo11N7m6+7UG^f1-n5(1>AvFl_Z18 z`{Q84^oDLg`rC`|cAU?v^}Nf%V{{Vi(k`D9JP@t}VCjv+GB4!aOhGRAaK><);?Nbo zDHPs@*94C-W$Z1$GT*(KritmRo6O=!b#g%w9#-3py>y{H4WsD?^nJl5$- zV|kkhtGx3HK7cfK^Q0ebUEPj%O%~eIc=YUj=H?G0L_`cM|6)rqWZvP?vpLTXF(J_T zN@5(T=s+_PVA&re)gnA$-jUAjc~!|)tkiTGt0LT@rO*XO7wl%zq6}u=d|DQ61jn`r zOG$_Ly3}FO4FwmFWLOPA?Dc+nZ>N8@?J0f59L{NN7y;8@f&wUwT;SHj9W&2!C)PVx zc_9EfRzTf{&eEWN8`G@oCWETD5+LcEeJ(_E@2*A#IZ*i>Z$z~6WgKD-MO^1<-mo29 z0Q7E(*{B)jkH{HPqFZ7q7vqXNxnvqV5t&zjLZ~|p zVS(;0`VMjJ1fR3;8{+TvmHUVlBuK#ga)VGebTu&{`-*jZ&LZ9mf#|w*pU9;GzQgkg z2pB3=%Te>)!wx8t>~e!)WP=hBYJ~p1GT_XmP2v7DFLbH>eu=X6b&R-biKGa*>|KrD zAVdx-Wp3@AG9uA@So*^tYUMpyk?LAuKZg-23&*k(N`qu_na`F@%@ z!SJ(Ws21nb^Ku ze6>UCj2e($OmdiS2% z9@efefUTu9?(O#FAqd%k<(a z>+nNfb=aq*Mo)53Pd78Bxvje66>zanldU8-M)(zA5DG<~`9js-(%t7O)oTk2OYdq6 ztX~_4mG|!JDT#ax}jtRY`BsMMhK=ypi+V1`*i^*cigxapHTn)DL_j1 z7YwJ~M@hX9q4&GU}Q?aAuSgukPLS^n;bu(g@MQf|-q2JW4aqL-Uy$wscd1 zp3(5_hAL{QrxAaM4!YU0PMLSmFRFK|8&H2hqtIy~^nU#`iy;UDrZ9hk*Zb*Gm~_tv zyDGDF@NOYNyb30|zlF4g6{_G#au?|8li@XBlsfj)Z?UgL3$Z5N=1}Vbe?9qpu#($4 zke9Iv;r&30g~KJilDB2`(_5GJ!iDJ=Dp!AKfLD*mg!ojk(vtxt;T@`9K5LOJR1ixt z94lA^D}!tVLIuULA073J?(GN=tS;#pwtq~Eb zUBtfYQ|{Z*X?go!lIlCF5r<3^%eV(ek+qEGX~F#rXh8ojMV0P(F+Vc&$W*sf5U(%E zKI1pRFLwNAz+mEip=+e>1W3V=?`d(S1-7^d{a4+kn4|y2&^?Qggb(8ig6FeP5Bz#Z zd>d|8Yx?2te~F0WHZz4x%P6=+Bc&>X#60Txj~uQ0IQE_;_77aSOsIFj_j*N|dQ;o% zcvkzFvG?|@*nVSvRo|O^3N)-l5QNyX5} ztGW_+w^Ez0eEQ;WE>r=jD8Dg?Ppf4gTb5D<&X{SgKh5av2JOw8;agqcey-f&1e9@U zepfJ8cUx|!6ZNnxsC@C*{*c3XvbVAQrg*tBc6eXkmlxUn-*{VD*DBO6^*=yK;o46; zPzg6K6!U@bSOlpXM;FR1)f;&iT7rz9ROVO3$TaUXOdCIo&Rqwbp8PO?Bf6?8Wer+D zxW-Z}=MKH6htNUCD~h+veM=hif?g#1AQf}1D+D=yKRqifHjm(Z>`jt((7*fZxpq$R zik2L27=E-Il?%dW;r4fNN9VYt0+p5oH9VSME%Ca|cAV$mrnV^lU9<1`16a?q8G|<| zyce_dlX+Ih^OS8~Op~nF1s6H1K$TZsv(l@V z3mXGpwNCvYU5uO=Z)0!MBWtd`)kZtk(dmPYnNtTF0>JL_jc`m0dk{e;i>oaWm%~Dr z$_-*)7bIt#NWS&0LHEvo`Fi>8%YV{cf3uyW zH&0`2iMGj0F{{AlsdVG&-9#`F_nHFt97hg)d=`)ORsoyCe@_x%zy(@@|0|fQ0R6t_ z?3o%FMZlPFe+0vr;3K{_3cV%@)FPF;fpxMskySEL_0hIgakB6uhYyBP-x^Z@$$&&a zD=Pz&Gks#B)MI0PV@$HiurR<%6rKNJYX3jXvHu?vUd{#H$;sG6*2me;-rnMaRxT89 yZiRDpV}bn(-y$#9ud$&G-r2dCjTx#h!z(LebDVs`pJu5Tm>GuQU?fd2C;tzt66c-( literal 0 HcmV?d00001 diff --git a/tests/golden/make_loudness_sqat.m b/tests/golden/make_loudness_sqat.m new file mode 100644 index 0000000..6517516 --- /dev/null +++ b/tests/golden/make_loudness_sqat.m @@ -0,0 +1,122 @@ +function make_loudness_sqat(sqatRoot) +% MAKE_LOUDNESS_SQAT Golden data for tests/test_psychoacoustics.py (loudness). +% +% make_loudness_sqat(sqatRoot) runs Loudness_ISO532_1 from SQAT +% (https://github.com/ggrecow/SQAT) on short synthetic signals and writes +% loudness_sqat.mat next to this file. Python rebuilds the same signals +% from the parameters stored here (see the helpers of the test file). +% +% Signals (sound pressure in Pa). Every signal is rounded to single +% precision before the call, as SignalObj stores float32, so MATLAB and +% Python see the same samples: +% tone a*sin(2*pi*f*t), a = sqrt(2)*2e-5*10^(L/20) +% multitone the 28 one-third-octave centre frequencies 25 Hz..12.5 kHz, +% each at L dB, Schroeder phases pi*k*(k-1)/28, summed in +% order k = 1..28 +% burst tone that is non-zero for t0 <= t < t1 only +% ramp tone whose level goes linearly from L to L2 dB over the +% signal: a(t) = sqrt(2)*2e-5*10.^((L+(L2-L)*t/T)/20) +% with t = (0:n-1)'/fs and T = n/fs. +% +% Stored per case (cell array "cases" of structs): the signal +% parameters, the call arguments (field 0 free / 1 diffuse, method 1 +% stationary / 2 time-varying, time_skip), sum and sum of squares of the +% single precision input (to check that Python rebuilt the same signal), +% and the outputs: stationary N, LN, Nspec [1x240]; time-varying N(t), LN(t), +% time, N5, Nmax, and the specific loudness at every 25th output frame +% (Nspec_t [frames x 240], frames in Nspec_idx, 1-based) to keep the file +% small. +% +% Also stored: matlab_version and sqat_commit (the commit of the SQAT +% checkout, read from its .git folder). +% +% Usage: matlab -singleCompThread -batch "make_loudness_sqat('/path/to/SQAT')" + +addpath(sqatRoot); +startup_SQAT([strip(sqatRoot, "right", filesep) filesep]); % expects a trailing separator + +defs = { +% name kind f L L2 t0 t1 fs dur field method skip + 'tone1k_40' 'tone' 1000 40 0 0 0 48000 1.0 0 1 0 + 'multitone45' 'multitone' 0 45 0 0 0 48000 1.0 0 1 0.2 + 'multitone45_D' 'multitone' 0 45 0 0 0 48000 1.0 1 1 0 + 'tone250_80_441' 'tone' 250 80 0 0 0 44100 1.0 0 1 0 + 'burst1k_70' 'burst' 1000 70 0 0.1 0.2 48000 0.5 0 2 0 + 'ramp1k_40_80_D' 'ramp' 1000 40 80 0 0 48000 1.0 1 2 0.2 + 'tone1k_60_441' 'tone' 1000 60 0 0 0 44100 0.5 0 2 0 +}; + +cases = cell(1, size(defs, 1)); % cell array: stationary and time-varying cases store different fields +for c = 1:size(defs, 1) + d = defs(c, :); + p = cell2struct(d, {'name','kind','f','L','L2','t0','t1','fs','dur', ... + 'field','method','time_skip'}, 2); + x = il_signal(p); + OUT = Loudness_ISO532_1(x, p.fs, p.field, p.method, p.time_skip, 0); + p.sum_x = sum(x); + p.sum_x2 = sum(x.^2); + if p.method == 1 + p.N = OUT.Loudness; + p.LN = OUT.LoudnessLevel; + p.Nspec = OUT.SpecificLoudness; + else + p.N = OUT.InstantaneousLoudness; + p.LN = OUT.InstantaneousLoudnessLevel; + p.time = OUT.time; + p.N5 = OUT.N5; + p.Nmax = OUT.Nmax; + idx = 1:25:numel(OUT.time); + p.Nspec_idx = idx; + p.Nspec_t = OUT.InstantaneousSpecificLoudness(idx, :); + end + cases{c} = p; + fprintf('%-16s N = %s\n', p.name, mat2str(p.N(1:min(3,end))', 6)); +end + +matlab_version = version; +sqat_commit = il_sqat_commit(sqatRoot); +out = fullfile(fileparts(mfilename('fullpath')), 'loudness_sqat.mat'); +save(out, 'cases', 'matlab_version', 'sqat_commit', '-v7'); +fprintf('wrote %s (%d cases)\n', out, numel(cases)); +end + +function x = il_signal(p) +n = round(p.dur * p.fs); +t = (0:n-1)' / p.fs; +switch p.kind + case 'tone' + x = sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*p.f*t); + case 'burst' + x = sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*p.f*t); + x(~(t >= p.t0 & t < p.t1)) = 0; + case 'ramp' + T = n / p.fs; + x = sqrt(2)*2e-5*10.^((p.L + (p.L2-p.L)*t/T)/20) .* sin(2*pi*p.f*t); + case 'multitone' + K = 28; + a = sqrt(2)*2e-5*10^(p.L/20); + x = zeros(n, 1); + for k = 1:K + fk = 1000 * 10^(((k-1)-16)/10); + x = x + a * sin(2*pi*fk*t + pi*k*(k-1)/K); + end +end +x = double(single(x)); +end + +function c = il_sqat_commit(sqatRoot) +% Commit of the SQAT checkout, read from its .git folder ('' if unknown). +c = ''; +g = fullfile(sqatRoot, '.git'); +if ~isfolder(g), return; end +head = strtrim(fileread(fullfile(g, 'HEAD'))); +if ~startsWith(head, 'ref: '), c = head; return; end +ref = head(6:end); +if isfile(fullfile(g, ref)) + c = strtrim(fileread(fullfile(g, ref))); +elseif isfile(fullfile(g, 'packed-refs')) + lines = splitlines(fileread(fullfile(g, 'packed-refs'))); + hit = lines(endsWith(lines, [' ' ref])); + if ~isempty(hit), c = extractBefore(hit{1}, ' '); end +end +end diff --git a/tests/golden/make_roughness_sqat.m b/tests/golden/make_roughness_sqat.m new file mode 100644 index 0000000..3549ce9 --- /dev/null +++ b/tests/golden/make_roughness_sqat.m @@ -0,0 +1,82 @@ +function make_roughness_sqat(sqatRoot) +% MAKE_ROUGHNESS_SQAT Golden data for tests/test_psychoacoustics.py (roughness). +% +% make_roughness_sqat(sqatRoot) runs Roughness_Daniel1997 from SQAT +% (https://github.com/ggrecow/SQAT) on short synthetic signals and writes +% roughness_sqat.mat next to this file. Python rebuilds the same signals +% from the parameters stored here (see the helpers of the test file). +% +% Signal (sound pressure in Pa), rounded to single precision before the +% call, as SignalObj stores float32: +% x = A*(1 + m*sin(2*pi*fm*t)).*sin(2*pi*fc*t), t = (0:n-1)'/fs, +% A = sqrt(2)*2e-5*10^(L/20)/sqrt(1 + m^2/2), +% so that L is the level of the rms of the modulated signal (m = 0 gives +% a pure tone). +% +% Stored per case (struct array "cases"): the signal parameters, +% time_skip, sum and sum of squares of the single precision input (to +% check that Python rebuilt the same signal), and the outputs R(t) +% (InstantaneousRoughness), time, the specific roughness +% (InstantaneousSpecificRoughness, 47 x frames) and Rmean. +% +% Also stored: matlab_version and sqat_commit (the commit of the SQAT +% checkout, read from its .git folder). +% +% Usage: matlab -singleCompThread -batch "make_roughness_sqat('/path/to/SQAT')" + +addpath(sqatRoot); +startup_SQAT([strip(sqatRoot, "right", filesep) filesep]); % expects a trailing separator + +defs = { +% name fc fm m L fs dur skip + 'am1k_70_ref' 1000 70 1 60 48000 2.0 0.5 + 'am1k_70_441' 1000 70 1 60 44100 1.0 0 + 'am2k_50_m05_32k' 2000 50 0.5 70 32000 1.0 0 + 'tone1k_60' 1000 0 0 60 48000 1.0 0 +}; + +cases = struct([]); +for c = 1:size(defs, 1) + p = cell2struct(defs(c, :), {'name','fc','fm','m','L','fs','dur','time_skip'}, 2); + n = round(p.dur * p.fs); + t = (0:n-1)' / p.fs; + A = sqrt(2)*2e-5*10^(p.L/20) / sqrt(1 + p.m^2/2); + x = double(single(A * (1 + p.m*sin(2*pi*p.fm*t)) .* sin(2*pi*p.fc*t))); + OUT = Roughness_Daniel1997(x, p.fs, p.time_skip, 0); + p.sum_x = sum(x); + p.sum_x2 = sum(x.^2); + p.R = OUT.InstantaneousRoughness; + p.time = OUT.time; + p.Rspec = OUT.InstantaneousSpecificRoughness; + p.Rmean = OUT.Rmean; + if isempty(cases) + cases = p; + else + cases(end+1) = p; %#ok + end + fprintf('%-16s Rmean = %.6f asper (%d frames)\n', p.name, p.Rmean, numel(p.R)); +end + +matlab_version = version; +sqat_commit = il_sqat_commit(sqatRoot); +out = fullfile(fileparts(mfilename('fullpath')), 'roughness_sqat.mat'); +save(out, 'cases', 'matlab_version', 'sqat_commit', '-v7'); +fprintf('wrote %s (%d cases)\n', out, numel(cases)); +end + +function c = il_sqat_commit(sqatRoot) +% Commit of the SQAT checkout, read from its .git folder ('' if unknown). +c = ''; +g = fullfile(sqatRoot, '.git'); +if ~isfolder(g), return; end +head = strtrim(fileread(fullfile(g, 'HEAD'))); +if ~startsWith(head, 'ref: '), c = head; return; end +ref = head(6:end); +if isfile(fullfile(g, ref)) + c = strtrim(fileread(fullfile(g, ref))); +elseif isfile(fullfile(g, 'packed-refs')) + lines = splitlines(fileread(fullfile(g, 'packed-refs'))); + hit = lines(endsWith(lines, [' ' ref])); + if ~isempty(hit), c = extractBefore(hit{1}, ' '); end +end +end diff --git a/tests/golden/roughness_sqat.mat b/tests/golden/roughness_sqat.mat new file mode 100644 index 0000000000000000000000000000000000000000..9b285f5eab1a734b9a9d9c0e89b76e1b7b8004c4 GIT binary patch literal 3725 zcma)L*ilHM)3C$yhz!ij`qDT{wCPk_S>489~ zN+;+gAeTcRQlbbaK!hWNkOx1!?}z8pvu5_p+H3aw)|x$Qula9Vk6r7h3C?{J)r?S9_MRGiJYpXkvw8s%}nl6UaSI?e8@ayiL> z6gQ4~8pG?k-ov)NmD6y8#b@HCmO%$zmz})w z=op8YY7s3J_5jrh!MAdAr{Q`utUpU75b5pIC}IJ%g*N%vOOSe;+lN?oV#TZPn5qC( zZR>^f;-FK5T>fL?kdQ?t%`Y|by`%j2c=;QSE1{RDib8V3>Hp*oww3$d?Y&LhWIWlX-lUqq(-D zg8TYVnXZ=i-_sV5Zdxr(c;@e?E40ntKs+Knv)xR+VL9i}j$ZGA|NTYpsnAlnR_{Ku zX^EDhrl!A{hIl@|9K_3KCz~A2lQ#@?T&k*;14U#V3$jVHpM7C7hp`N(|{ky!WmpW>%_B8vb_Mu9w59PBMPrX`qY>IbT>BD$lRXt7$VtlTQIxFh@Ji;%{|%4e%G9&-V;u-UY`Q_Xl4IA11X(N*i7 zW>U4dhyB?a8l_r4Vg1K=W_qf?do*IBJ&x zJM@%dHwSUU4{%yt%mbeoVmQX#C#M?|B4mKH>jLt)4vS}i7}$OuP{D~OqC^%A+2GA! zbi{G2?tq63)P?2&pHd1=z;BT-27eIBLek1W6a{Pa4Je7U{IH_|_E_HgN89{!VujsyyAqT13+ z(m-<)cG&w)fhwNSgXY&Dcf_MqK$_#vs;<@5%grke`=y9d3XI=Z%i*dkqDTq<62U~q zVGuJ%_(1$wpQePD1=6$yLJi2)73k;u*@gxD&AGowlQn(9-(eA+V1^J^_;g zT#g48J^hm;`y5WZY{KPKV;r1>KczwbK~-Nu`cjcDF% zpa0u?M}yALJ?_I`1Z53QZy5L!nwdx0*jTKUDJ~mlSHYue)A7f&dFP)@uQgUCa}N4& zm*2NXOP+8%&7ad_Rr1gElc&w{lCk{4NBVt=PnFm?3ekDR=aSeDh@U&x6X1D5Eb4If z<12qDJ{g^vX2|~tDgE@b+hoq#ZnWLDBae)*BNV7irhF^rQXkHdR1*w(<>b9=f9gq{ zH#Yy3R}w9!UCw^^X0%iHhS@>7{^5mMqSpH5={-Bn=n}F7@nB)S!4V zW%ZdOJnHG8drxA{t6!5Dzbfazb5&+3`Px#}t86-yh{acVeQGJ=6Uz_q3-GE}O66{s z!><%NT21`^*Tz1g-KwFlBqFzI&uS|u=&GBZnT1vVOUe$kiwa-gFeyBgE zPyawT*e} zAbjUmeq{m!ZrdTt{h61;Xo z3s)RRc=q>6t>4y={;bH!5S~S z-xg7lvoX;fubff(OEj9g@h+jX)OxQ{TzHFbxi7N{7*w-FP!FbdL<;p5oooFN6D!OS z*Ba^WVW+P;)luu2f88?1&)<*&l=Oc#q+fG=i%AQs&&p8vw&1(;B8>5|)iUy9(dV~C z)ZicBg23~*>5r)0*%QHHn!aI4j8AQqpvC}eg>xs3^N!l4vp3`zwzd14W4L?Vs>+iJPT-=V%2{bqm1Xvwuu zEVPc-%fmnCTpJpGRt|I;2tt;~=0rw0vhV*~m7Y(`cnL zkEpf2VM-ojbPq)UYimgHThQVHcaPcn&@%*7Vd9I+R&+Bd-OWNO0mZ4th?Tfs(R170 zRVA6DTNtaK(rz!G3G6hoz5^`lIQ@R}T+j#Zh5XaE4()ciVLit}djKE3w1y7DJl`1OD=sxygujx^ zTb+9$ReQ&yQMEKdMXP*MsoylXd78*_Xg6%rS{Bl!+bc`PIT1J z7S|#ATV)_E276P>q3L92Z~{XHRHwvpKq9rXv>p!uyBYCGXE*{`*z64S*o)+fjxFu0 znc|q)eDHSl7bnXZ28D16I*yz@3=Hru@9{w#SC$QZiy@(s@pkR8)q7VybS)~Uw++AK z@u@rOSKQO`dS3!zB3x3Lh8(aA?r!a@eC_N=V{|#siY9(vd??#`enIhrE|`cY6>x|G zaljblq?rRyb7VOLgAygs?45}sK0*8nqcc;i($}Lm5RX`|XLqWoTz~#mHv8*--n@^e zXEK|&rosW#ZSb z-NR#@CdShJ{hj@amKg`wZdY@;|6#iQKePY8nTM<*4uyqzd0WPW$A*OXb6n5g2Y&f_ x=`%R4rm92I))?sT1+_na83V_pR6l= p.t0) & (t < p.t1))] = 0 + elif p.kind == 'ramp': + T = n / p.fs + x = (np.sqrt(2) * 2e-5 * 10 ** ((p.L + (p.L2 - p.L) * t / T) / 20) + * np.sin(2 * np.pi * p.f * t)) + elif p.kind == 'multitone': + K = 28 + a = np.sqrt(2) * 2e-5 * 10 ** (p.L / 20) + x = np.zeros(n) + for k in range(1, K + 1): + fk = 1000 * 10 ** (((k - 1) - 16) / 10) + x = x + a * np.sin(2 * np.pi * fk * t + np.pi * k * (k - 1) / K) + return x.astype(np.float32) + + +def roughness_signal(p): + """Rebuild the signal of one case of make_roughness_sqat.m.""" + n = int(round(p.dur * p.fs)) + t = np.arange(n) / p.fs + A = np.sqrt(2) * 2e-5 * 10 ** (p.L / 20) / np.sqrt(1 + p.m ** 2 / 2) + x = A * (1 + p.m * np.sin(2 * np.pi * p.fm * t)) * np.sin(2 * np.pi * p.fc * t) + return x.astype(np.float32) + + +def same_samples(tc, p, x): + """The rebuilt float32 signal is the one MATLAB used.""" + x = x.astype(float) + tc.assertLess(abs(np.sum(x ** 2) / p.sum_x2 - 1), 1e-12) + tc.assertLess(abs(np.sum(x) - p.sum_x), 1e-12 * np.sum(np.abs(x))) + + +def band_noise(seed, f1, f2, level, fs=48000, dur=2.0): + """Noise with a flat spectrum between f1 and f2 (Hz) at level dB SPL.""" + n = int(fs * dur) + f = np.fft.rfftfreq(n, 1 / fs) + rng = np.random.default_rng(seed) + spec = np.exp(2j * np.pi * rng.random(f.size)) * ((f >= f1) & (f <= f2)) + x = np.fft.irfft(spec, n) + return (x * P0 * 10 ** (level / 20) / np.sqrt(np.mean(x ** 2))).astype(np.float32) + + +def at_loudness(target, x, fs=48000): + """Scale x so that its stationary free-field loudness is target sone.""" + for _ in range(3): + sig = pytta.SignalObj(x, 'time', fs) + N = pa.loudness_zwicker(sig)['N'] + x = (x * 10 ** (10 * np.log2(target / N) / 20)).astype(np.float32) + return x + + +class TestLoudness(unittest.TestCase): + + def test_1khz_tone_at_40_db_is_1_sone_and_40_phon(self): + # Definition of the sone (ISO 532-1, 3.13 and Annex C). Zwicker's + # tables reproduce it to 2.4e-5 sone (SQAT gives 0.999976). + fs = 48000 + t = np.arange(fs) / fs + sig = pytta.SignalObj(np.sqrt(2) * P0 * 100 * np.sin(2 * np.pi * 1000 * t), + 'time', fs) + out = pa.loudness_zwicker(sig) + print(f'\n N = {out["N"]:.6f} sone, LN = {out["LN"]:.4f} phon') + self.assertAlmostEqual(out['N'], 1, delta=1e-4) + self.assertAlmostEqual(out['LN'], 40, delta=1e-3) + self.assertEqual(out['specific'].shape, (240,)) + np.testing.assert_allclose(out['bark'], np.arange(1, 241) / 10) + + def test_annex_b_stationary_reference_values(self): + # ISO 532-1:2017 Annex B.2 and B.3, values of the test implementation + # (test signals 1 to 4). Conformance (5.1): within 5 % or 0.1 sone. + # Signal 1 is given as one-third-octave levels; 2 to 4 are tones, + # synthesised here at the nominal level and analysed from 0.2 s on, + # as B.1 prescribes (the filter transients at the start add up to + # 1.6 % otherwise). + levels = [-60, -60, 78, 79, 89, 72, 80, 89, 75, 87, 85, 79, 86, 80, + 71, 70, 72, 71, 72, 74, 69, 65, 67, 77, 68, 58, 45, 30] + out = pa._loudness_from_third_octave_levels(levels, 'free') + N, LN = out['N'], out['LN'] + print(f'\n signal 1: N = {N:.3f} (ref 83.296), LN = {LN:.3f} (ref 103.802)') + self.assertLessEqual(abs(N - 83.296), max(0.05 * 83.296, 0.1)) + fs = 48000 + t = np.arange(fs) / fs + for f, L, ref in [(250, 80, 14.655), (1000, 60, 4.019), (4000, 40, 1.549)]: + sig = pytta.SignalObj(np.sqrt(2) * P0 * 10 ** (L / 20) + * np.sin(2 * np.pi * f * t), 'time', fs) + N = pa.loudness_zwicker(sig, time_skip=0.2)['N'] + print(f' {f} Hz {L} dB: N = {N:.4f} sone, ref {ref}, diff {N - ref:+.4f}') + self.assertLessEqual(abs(N - ref), max(0.05 * ref, 0.1)) + + def test_matches_sqat_stationary(self): + # Same float32 samples in both; what is left is double rounding in + # the filters (and in the 44.1 to 48 kHz resampling, which follows + # MATLAB's resample to 2e-15), far below 1e-9 sone. + for p in golden_cases('loudness_sqat.mat'): + if p.method != 1: + continue + x = loudness_signal(p) + same_samples(self, p, x) + sig = pytta.SignalObj(x, 'time', p.fs) + out = pa.loudness_zwicker(sig, field=['free', 'diffuse'][p.field], + time_skip=p.time_skip) + dN = abs(out['N'] - p.N) + dLN = abs(out['LN'] - p.LN) + dspec = np.max(np.abs(out['specific'] - p.Nspec)) + print(f'\n {p.name}: N = {out["N"]:.6f}, |dN| = {dN:.1e} sone, ' + f'|dLN| = {dLN:.1e} phon, max |dN\'| = {dspec:.1e} sone/Bark') + self.assertLess(dN, 1e-9) + self.assertLess(dLN, 1e-9) + self.assertLess(dspec, 1e-9) + + def test_matches_sqat_time_varying(self): + # The nonlinear decay of step 8 picks its branch from comparisons + # between successive values: core loudness that differ by 3e-14 + # sone (the rounding of the filter bank, MATLAB against scipy) come + # out of it up to 5e-7 sone apart on the decay after a tone burst. + # Tolerance 1e-6 sone (1e-5 phon). + for p in golden_cases('loudness_sqat.mat'): + if p.method != 2: + continue + x = loudness_signal(p) + same_samples(self, p, x) + sig = pytta.SignalObj(x, 'time', p.fs) + out = pa.loudness_zwicker(sig, field=['free', 'diffuse'][p.field], + method='time_varying', + time_skip=p.time_skip) + self.assertEqual(out['N'].shape, p.N.shape) + self.assertEqual(out['specific'].shape, (p.N.size, 240)) + np.testing.assert_allclose(out['time'], p.time, rtol=0, atol=1e-12) + dN = np.max(np.abs(out['N'] - p.N)) + dLN = np.max(np.abs(out['LN'] - p.LN)) + dspec = np.max(np.abs(out['specific'][p.Nspec_idx - 1] - p.Nspec_t)) + dstat = max(abs(out['N5'] - p.N5), abs(out['Nmax'] - p.Nmax)) + print(f'\n {p.name}: Nmax = {out["Nmax"]:.5f}, max |dN(t)| = {dN:.1e}, ' + f'|dLN| = {dLN:.1e}, |dN\'| = {dspec:.1e}, |dN5|,|dNmax| <= {dstat:.1e}') + self.assertLess(dN, 1e-6) + self.assertLess(dLN, 1e-5) + self.assertLess(dspec, 1e-6) + self.assertLess(dstat, 1e-6) + + def test_rejects_invalid_input(self): + fs = 48000 + two = pytta.SignalObj(np.zeros((fs, 2)), 'time', fs) + one = pytta.SignalObj(np.zeros(fs), 'time', fs) + with self.assertRaises(ValueError): + pa.loudness_zwicker(two) + with self.assertRaises(ValueError): + pa.loudness_zwicker(one, field='near') + with self.assertRaises(ValueError): + pa.loudness_zwicker(one, method='fast') + with self.assertRaises(ValueError): + pa.loudness_zwicker(one, time_skip=1.0) # not shorter than the signal + + +class TestSharpness(unittest.TestCase): + + def test_din_reference_noise_is_1_acum(self): + # DIN 45692, 3.1 and 6: narrowband noise 920 Hz to 1080 Hz at 60 dB + # gives 1.00 acum. Tolerance: the DIN limit for measuring systems, + # 5 % or 0.05 acum. This noise (seed 7) is not the one k was fitted on. + x = band_noise(7, 920, 1080, 60) + sig = pytta.SignalObj(x, 'time', 48000) + S = pa.sharpness_din45692(sig)['S'] + print(f'\n S = {S:.4f} acum, diff {S - 1:+.4f}') + self.assertLessEqual(abs(S - 1), 0.05) + + def test_din_table_a2_a3_values(self): + # DIN 45692 Tables A.2 (critical-band noises) and A.3 (broadband noise + # from fu to 10 kHz), all at the loudness of the reference signal, + # 4 sone. Tolerance: 5 % or 0.05 acum (DIN 45692, 6). + cases = [ # (f1, f2, reference acum) + (200, 300, 0.38), (920, 1080, 1.00), (3700, 4400, 2.82), + (9500, 12000, 8.55), + (250, 10000, 2.70), (1000, 10000, 3.05), (4000, 10000, 4.49), + (8500, 10000, 7.46)] + for f1, f2, ref in cases: + x = at_loudness(4.0, band_noise(11, f1, f2, 60)) + sig = pytta.SignalObj(x, 'time', 48000) + S = pa.sharpness_din45692(sig)['S'] + print(f'\n {f1}-{f2} Hz: S = {S:.3f} acum, ref {ref}, diff {S - ref:+.3f}', + end='') + self.assertLessEqual(abs(S - ref), max(0.05 * ref, 0.05)) + + def test_time_varying_s50_of_reference_noise_is_1_acum(self): + # DIN 45692, 6, note: S50 is the value to check with narrowband + # noise, whose sharpness fluctuates. Same tolerance as above. + x = band_noise(7, 920, 1080, 60, dur=1.0) + sig = pytta.SignalObj(x, 'time', 48000) + out = pa.sharpness_din45692(sig, method='time_varying', time_skip=0.2) + print(f'\n S50 = {out["S50"]:.4f} acum over {out["S"].size} frames') + self.assertEqual(out['S'].shape, out['time'].shape) + self.assertLessEqual(abs(out['S50'] - 1), 0.05) + + def test_silence_has_zero_sharpness(self): + # N' = 0 everywhere leaves Formula (1) as 0/0; it is reported as 0. + sig = pytta.SignalObj(np.zeros(48000), 'time', 48000) + self.assertEqual(pa.sharpness_din45692(sig)['S'], 0) + + def test_rejects_two_channels(self): + sig = pytta.SignalObj(np.zeros((48000, 2)), 'time', 48000) + with self.assertRaises(ValueError): + pa.sharpness_din45692(sig) + + +class TestRoughness(unittest.TestCase): + + def test_reference_am_tone_is_1_asper(self): + # Definition of the asper: 1 kHz tone, 100 % amplitude modulated at + # 70 Hz, 60 dB (rms of the modulated signal). SQAT calibrates on its + # reference .wav file with a 7-digit constant and gives 1.000001 for + # this synthetic version; 1e-4 asper leaves room for that only. + fs = 48000 + t = np.arange(2 * fs) / fs + A = np.sqrt(2) * P0 * 1000 / np.sqrt(1.5) + x = A * (1 + np.sin(2 * np.pi * 70 * t)) * np.sin(2 * np.pi * 1000 * t) + out = pa.roughness_daniel1997(pytta.SignalObj(x, 'time', fs), time_skip=0.5) + print(f'\n Rmean = {out["Rmean"]:.6f} asper, diff {out["Rmean"] - 1:+.1e}') + self.assertAlmostEqual(out['Rmean'], 1, delta=1e-4) + self.assertEqual(out['specific'].shape, (out['R'].size, 47)) + np.testing.assert_allclose(out['bark'], np.arange(1, 48) / 2) + + def test_matches_sqat(self): + # The modulation depth and the correlation of channels whose + # envelope sits at the FFT rounding floor depend on the rounding + # pattern of the FFT (SQAT notes up to 5e-4 asper between two FFT + # orderings in MATLAB); numpy's FFT is another ordering. Tolerance + # 1e-3 asper on R(t) and on the specific roughness. + for p in golden_cases('roughness_sqat.mat'): + x = roughness_signal(p) + same_samples(self, p, x) + sig = pytta.SignalObj(x, 'time', p.fs) + out = pa.roughness_daniel1997(sig, time_skip=p.time_skip) + self.assertEqual(out['R'].shape, p.R.shape) + np.testing.assert_allclose(out['time'], p.time, rtol=0, atol=1e-12) + dR = np.max(np.abs(out['R'] - p.R)) + dspec = np.max(np.abs(out['specific'] - p.Rspec.T)) + dmean = abs(out['Rmean'] - p.Rmean) + print(f'\n {p.name}: Rmean = {out["Rmean"]:.6f}, max |dR(t)| = {dR:.1e}, ' + f'max |dR\'| = {dspec:.1e}, |dRmean| = {dmean:.1e} asper') + self.assertLess(dR, 1e-3) + self.assertLess(dspec, 1e-3) + self.assertLess(dmean, 1e-3) + + def test_rejects_invalid_input(self): + fs = 48000 + two = pytta.SignalObj(np.zeros((fs, 2)), 'time', fs) + with self.assertRaises(ValueError): + pa.roughness_daniel1997(two) + short = pytta.SignalObj(np.zeros(fs // 10), 'time', fs) + with self.assertRaises(ValueError): + pa.roughness_daniel1997(short) # shorter than one 200 ms window + + +if __name__ == '__main__': + unittest.main()