From 1df4465bcda63223ab24b3a455735bba3dc4d347 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 21:39:39 +0100 Subject: [PATCH 1/4] [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() From a9224ddfcd56303a50288c17fddc1c290e3a2c85 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 23:40:37 +0100 Subject: [PATCH 2/4] [FEATURE] add the tonality of Aures to pytta.psychoacoustics tonality_aures1985(signal, field='free', time_skip=0.0) gives the tonality K(t) and its mean in tonality units: the tonal components of Terhardt, Stoll and Seewann (1982) in Hann frames of 160 ms, weighted by bandwidth, frequency and level excess, and the loudness ratio WGr of the signal without and with them (Aures 1985). It is written from the two papers. SQAT's Tonality_Aures1985 is GPL-3.0 and served only as the test oracle. The private helper of the Terhardt excitation is shared with the roughness port, whose parity with SQAT is unchanged. Tests check the unit (a 1 kHz tone at 60 dB is 1 t.u.), the level excess weighting against its formula, parity with SQAT on golden data from tests/golden/make_tonality_sqat.m (tones within 3.4e-5 t.u., a tone in noise within 0.012 t.u.) and the rejection of invalid input. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/index.rst | 5 +- pytta/__init__.py | 5 +- pytta/psychoacoustics.py | 217 ++++++++++++++++++++++++++++-- tests/golden/make_tonality_sqat.m | 100 ++++++++++++++ tests/golden/tonality_sqat.mat | Bin 0 -> 1053 bytes tests/test_psychoacoustics.py | 79 ++++++++++- 6 files changed, 392 insertions(+), 14 deletions(-) create mode 100644 tests/golden/make_tonality_sqat.m create mode 100644 tests/golden/tonality_sqat.mat diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ab5f549..acf770c 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -58,8 +58,9 @@ Modules: 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; + loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel + and Weber 1997) and tonality (Aures 1985) of a sound pressure + signal; .. toctree:: classes/index diff --git a/pytta/__init__.py b/pytta/__init__.py index a483f56..6f4b03e 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -73,8 +73,9 @@ 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; + loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel + and Weber 1997) and tonality (Aures 1985) of a sound pressure + signal; You can find usage examples in the toolbox's examples folder. diff --git a/pytta/psychoacoustics.py b/pytta/psychoacoustics.py index 58c4e01..8f12817 100644 --- a/pytta/psychoacoustics.py +++ b/pytta/psychoacoustics.py @@ -9,11 +9,14 @@ >>> pytta.psychoacoustics.loudness_zwicker(signal) >>> pytta.psychoacoustics.sharpness_din45692(signal) >>> pytta.psychoacoustics.roughness_daniel1997(signal) + >>> pytta.psychoacoustics.tonality_aures1985(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). +* roughness_daniel1997: roughness after Daniel and Weber (1997); +* tonality_aures1985: tonality after Aures (1985), tonal components after + Terhardt, Stoll and Seewann (1982). 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 @@ -53,6 +56,11 @@ sharpness_din45692 is written from the text of DIN 45692:2009 (Formula (1)). +tonality_aures1985 is written from the papers cited below. SQAT's +Tonality_Aures1985 is GPL-3.0; it was run only as a test oracle, and of its +file only the help text (calling convention) was read. Where the papers +leave a choice open, the function docstring says what was chosen. + References ---------- @@ -66,6 +74,16 @@ Daniel, P., & Weber, R. (1997). Psychoacoustical roughness: implementation of an optimized model. Acustica, 83, 113-123. +Aures, W. (1985). Berechnungsverfahren fuer den sensorischen Wohlklang +beliebiger Schallsignale. Acustica, 59, 130-141. + +Terhardt, E., Stoll, G., & Seewann, M. (1982). Algorithm for extraction of +pitch and pitch salience from complex tonal signals. The Journal of the +Acoustical Society of America, 71(3), 679-688. doi:10.1121/1.387544 + +Fastl, H., & Zwicker, E. (2007). Psychoacoustics: Facts and Models (3rd ed.). +Springer. + 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, @@ -139,7 +157,8 @@ import numpy as np import scipy.signal as ss -__all__ = ['loudness_zwicker', 'sharpness_din45692', 'roughness_daniel1997'] +__all__ = ['loudness_zwicker', 'sharpness_din45692', 'roughness_daniel1997', + 'tonality_aures1985'] _FIELDS = ('free', 'diffuse') _METHODS = ('stationary', 'time_varying') @@ -777,8 +796,8 @@ def _mround(v): 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.""" +def _excitation_params(N, fs): + """Bins, Bark scale, thresholds and a0 of the Terhardt excitation 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 @@ -795,6 +814,13 @@ def _roughness_params(N, fs): 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) + return {'N': N, 'n01': n01, 'bins': bins, 'freqs': freqs, 'bark': bark, + 'min_exc': min_exc, 'min_bf': min_bf, 'a0': a0} + + +def _roughness_params(N, fs): + """_excitation_params plus H(fmod) and g(z) of the roughness model.""" + prm = _excitation_params(N, fs) hw = np.zeros((_CHNO, N)) for proto, fmax, rows in ((_H2, 358, (2,)), (_H5, 502, (5,)), (_H16, 645, (16,)), (_H21, 645, (21,)), @@ -808,10 +834,9 @@ def _roughness_params(N, fs): 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} + prm['hw'] = hw + prm['gzi'] = np.sqrt(np.interp(np.arange(1, _CHNO + 1) / 2, np.arange(25), _GR)) + return prm def _terhardt_excitation(spec, prm): @@ -965,3 +990,179 @@ def roughness_daniel1997(signal, time_skip=0.0): return {'R': R, 'specific': spec, 'time': time, 'bark': np.arange(1, _CHNO + 1) / 2, 'Rmean': float(np.mean(R[_skip_index(time, time_skip):]))} + + +# --------------------------------------------------------------------------- +# Tonality, Aures (1985), tonal components after Terhardt et al. (1982) +# --------------------------------------------------------------------------- + +_T_WINDOW = 0.16 # s, Hann frames, 50 % overlap +_T_FMIN, _T_FMAX = 20.0, 5000.0 # Hz, range of the tonal components +_T_NEIGHBOURS = (25.0, 37.5) # Hz, Terhardt's samples i +- 2, i +- 3 at 12.5 Hz +_T_JUMP = 7.0 # dB, Terhardt et al. (1982), Eq. (2) +# c of Aures' Eq. (13), set so that a 1 kHz tone at 60 dB gives 1 t.u. (Aures +# 3.1); about 1.09 in an ideal implementation (Hastings et al. 2003) +_C_K = 1.105464 + + +def _bark_terhardt(f): + """Critical band rate (Bark) of f (Hz), Terhardt et al. (1982), Eq. (6).""" + f = np.asarray(f, dtype=float) / 1000 + return 13 * np.arctan(0.76 * f) + 3.5 * np.arctan((f / 7.5) ** 2) + + +def _threshold_terhardt(f): + """Threshold in quiet (dB SPL) of f (Hz), Terhardt et al. (1982), Eq. (8).""" + f = np.asarray(f, dtype=float) / 1000 + return 3.64 * f ** -0.8 - 6.5 * np.exp(-0.6 * (f - 3.3) ** 2) + 1e-3 * f ** 4 + + +def _tonal_components(P, df): + """ + Tonal components of one power spectrum P (Pa^2 per bin, summing to the + power of a component over its main lobe). Returns the peak bins, the + frequencies (Hz), levels (dB SPL) and SPL excess (dB) of every + component found, and the half-width (bins) of the group of samples + that represents a component. + """ + L = 10 * np.log10(P / _I_REF + _TINY) + j2, j3 = (_mround(h / df) for h in _T_NEIGHBOURS) + i = np.arange(max(j3 + 1, _mround(_T_FMIN / df)), + min(_mround(_T_FMAX / df), P.size - 2 - j3) + 1) + peak = (L[i - 1] < L[i]) & (L[i] >= L[i + 1]) + for j in (-j3, -j2, j2, j3): + peak &= L[i] - L[i + j] >= _T_JUMP + i = i[peak] + if i.size == 0: + return i, np.zeros(0), np.zeros(0), np.zeros(0), j3 + # frequency: Terhardt's Eq. (3), 0.46 Hz/dB at his 12.5 Hz spacing + fc = (i + 0.46 / 12.5 * (L[i + 1] - L[i - 1])) * df + # level: power of the main lobe of the Hann window (5 samples) + lc = 10 * np.log10(sum(P[i + d] for d in range(-2, 3)) / _I_REF) + # SPL excess, Terhardt et al. (1982), Eqs. (4), (5) and (7); the noise + # skips the central samples of every component (25 Hz on each side) + zc = _bark_terhardt(fc) + zb = _bark_terhardt(np.arange(P.size) * df) + noisy = np.ones(P.size, dtype=bool) + for d in range(-j2, j2 + 1): + noisy[i + d] = False + lx = np.zeros(i.size) + for u in range(i.size): + lower = fc[u] <= fc # f_u below f_v: lower slope of v + s = np.where(lower, 27.0, + np.minimum(-24 - 230 / fc + 0.2 * lc, 0.0)) + le = lc - s * (zc - zc[u]) + others = np.arange(i.size) != u + tones = np.sum(10 ** (le[others] / 20)) ** 2 + noise = np.sum(P[noisy & (np.abs(zb - zc[u]) <= 0.5)]) / _I_REF + lx[u] = lc[u] - 10 * np.log10(tones + noise + + 10 ** (_threshold_terhardt(fc[u]) / 10)) + return i, fc, lc, lx, j3 + + +def _third_octave_levels(P, df): + """Levels (dB SPL) of the 28 one-third-octave bands of a power spectrum.""" + f = np.arange(P.size) * df + out = np.empty(28) + for n, fm in enumerate(_OB_FC): + band = (f >= fm * 10 ** -0.05) & (f < fm * 10 ** 0.05) + out[n] = 10 * np.log10(np.sum(P[band]) / _I_REF + _TINY) + return out + + +def tonality_aures1985(signal, field='free', time_skip=0.0): + """ + Tonality after Aures (1985), with the tonal components of Terhardt, + Stoll and Seewann (1982). + + Written from the two papers; SQAT's Tonality_Aures1985 is GPL-3.0 and + served only as a test oracle. Per Hann frame of 160 ms (50 % overlap): + + * tonal components (Terhardt et al., Eqs. 1 to 3): local maxima from + 20 Hz to 5 kHz at least 7 dB above the samples 25 and 37.5 Hz away + on both sides (Terhardt's i +- 2 and i +- 3 at his 12.5 Hz spacing), + frequency by Eq. (3), level from the power of the main lobe; + * SPL excess LX (Eqs. 4 to 8): level minus the masking by the other + components (triangular excitation, slopes 27 dB/Bark and + 24 + 230/f - 0.2 L dB/Bark), by the noise of the critical band around + the component (samples within 25 Hz of a component left out) and by + the threshold in quiet; only components with LX > 0 count; + * WT = sqrt(sum (w2' w3')^2) (Aures, Eqs. 9 to 11), w2' the frequency + weighting (maximum at 700 Hz), w3' = 1 - exp(-LX / 15 dB); + * WGr = 1 - N_Gr / N (Eq. 12): ISO 532-1 loudness (stationary, from the + one-third-octave levels of the spectrum) of the spectrum without the + relevant components (37.5 Hz on each side replaced by the mean of + the next samples) over that of the whole spectrum; + * K = c WT^0.29 WGr^0.79 (Eq. 13), c = 1.105 so that a 1 kHz tone at + 60 dB gives 1. + + The components found this way are sinusoids, so the bandwidth + weighting w1 of Aures (Eq. 7) is 1. Aures also extracts narrowband + noise regions (narrower than a critical band, 7 dB above the + neighbouring bands), the only components with a bandwidth; that step + is not described precisely enough to reproduce and is left out, so + narrowband noise is rated by the peaks of its spectrum. Against SQAT: + within 1e-4 for tones, 0.015 for a tone in broadband noise. + + Parameters (default), (type): + ----------------------------- + + * signal (), (SignalObj): + One channel, sound pressure in Pa, at least 160 ms long, + sampled above 10 kHz; + + * field ('free'), (str): + 'free' or 'diffuse', for the loudness of WGr; + + * time_skip (0.0), (float): + Seconds excluded from Kmean at the start. + + Returns: + -------- + + dict, nt frames every 80 ms: + + * 'K': tonality K(t), t.u., shape (nt,); + * 'WT', 'WGr': tonal and loudness weightings, shape (nt,); + * 'time': start of each frame, s, shape (nt,); + * 'Kmean': mean of K(t) from the frame closest to time_skip on, + t.u. (float). + + """ + if field not in _FIELDS: + raise ValueError(f"field must be one of {_FIELDS}, got {field!r}.") + x, fs = _mono(signal) + if fs <= 2 * _T_FMAX: + raise ValueError(f'the sampling rate ({fs} Hz) must be above ' + f'{2 * _T_FMAX:.0f} Hz.') + N = 2 * _mround(fs * _T_WINDOW / 2) + hop = N // 2 + if x.size < N: + raise ValueError(f'the signal ({x.size / fs:.3f} s) is shorter than ' + f'one frame ({_T_WINDOW} s).') + nwin = (x.size - N) // hop + 1 + w = ss.windows.hann(N, sym=False) + scale = 2 / (np.sum(w) ** 2 * (N * np.sum(w ** 2) / np.sum(w) ** 2)) + df = fs / N + WT = np.zeros(nwin) + WGr = np.zeros(nwin) + for n in range(nwin): + X = np.fft.rfft(x[n * hop:n * hop + N] * w) + P = scale * np.abs(X) ** 2 + P[0] /= 2 + P[-1] /= 2 + i, fc, lc, lx, j3 = _tonal_components(P, df) + rel = lx > 0 + w2 = 1 / np.sqrt(1 + 0.2 * (fc[rel] / 700 + 700 / fc[rel]) ** 2) + w3 = 1 - np.exp(-lx[rel] / 15) + WT[n] = np.sqrt(np.sum((w2 * w3) ** 2)) + noise = P.copy() + for k in i[rel]: + noise[k - j3:k + j3 + 1] = (P[k - j3 - 1] + P[k + j3 + 1]) / 2 + Nt = _loudness_from_third_octave_levels(_third_octave_levels(P, df), field)['N'] + Ng = _loudness_from_third_octave_levels(_third_octave_levels(noise, df), field)['N'] + WGr[n] = max(1 - Ng / Nt, 0.0) if Nt > 0 else 0.0 + K = _C_K * WT ** 0.29 * WGr ** 0.79 + time = np.arange(nwin) * hop / fs + return {'K': K, 'WT': WT, 'WGr': WGr, 'time': time, + 'Kmean': float(np.mean(K[_skip_index(time, time_skip):]))} diff --git a/tests/golden/make_tonality_sqat.m b/tests/golden/make_tonality_sqat.m new file mode 100644 index 0000000..b57ec0a --- /dev/null +++ b/tests/golden/make_tonality_sqat.m @@ -0,0 +1,100 @@ +function make_tonality_sqat(sqatRoot) +% MAKE_TONALITY_SQAT Golden data for tests/test_psychoacoustics.py (tonality). +% +% make_tonality_sqat(sqatRoot) runs Tonality_Aures1985 from SQAT +% (https://github.com/ggrecow/SQAT) on short synthetic signals and writes +% tonality_sqat.mat next to this file. SQAT is used here only as a test +% oracle: that SQAT file is GPL-3.0 and pytta.psychoacoustics is written +% from the papers. Python rebuilds the same signals from the +% parameters stored here (see the helpers of the test file). +% +% Signals (sound pressure in Pa), rounded to single precision before the +% call, as SignalObj stores float32, t = (0:n-1)'/fs: +% 'tone': sum over the frequencies f of sqrt(2)*2e-5*10^(L/20)*sin(2*pi*f*t); +% 'tone+noise': the tone plus a multisine noise with components every +% 1 Hz from f1 to f2, phases 2*pi*mod(g*k^2, 1) with +% g = 0.6180339887 (k = 1..K), which scatter the group delay +% over the 1 s period so that every frame holds the whole band, +% scaled to the rms level Ln. +% +% Tonality: SQAT with the free-field loudness, time_skip 0. +% Stored per case (struct array "cases"): the signal parameters, the sum +% and sum of squares of the single precision input, and the outputs K(t), +% time and Kmean. Also stored: matlab_version and sqat_commit. +% +% Usage: matlab -singleCompThread -batch "make_tonality_sqat('/path/to/SQAT')" + +addpath(sqatRoot); +startup_SQAT([strip(sqatRoot, "right", filesep) filesep]); % expects a trailing separator + +tdefs = { +% name kind f L f1 f2 Ln fs dur + 'tone1k_60' 'tone' 1000 60 0 0 0 48000 1.0 + 'tone500_60' 'tone' 500 60 0 0 0 48000 1.0 + 'tone4k_60' 'tone' 4000 60 0 0 0 48000 1.0 + 'tone1k_30' 'tone' 1000 30 0 0 0 48000 1.0 + 'tones1k_1k2' 'tone' [1000 1200] 60 0 0 0 48000 1.0 + 'tone1k_noise60' 'tone+noise' 1000 60 200 5000 60 48000 1.0 + 'tone1k_noise50' 'tone+noise' 1000 60 200 5000 50 48000 1.0 +}; +cases = struct([]); +for c = 1:size(tdefs, 1) + p = cell2struct(tdefs(c, :), {'name','kind','f','L','f1','f2','Ln','fs','dur'}, 2); + t = (0:round(p.dur * p.fs)-1)' / p.fs; + x = zeros(size(t)); + for f = p.f + x = x + sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*f*t); + end + if strcmp(p.kind, 'tone+noise') + K = p.f2 - p.f1 + 1; + nz = zeros(size(t)); + for k = 1:K + nz = nz + sin(2*pi*(p.f1 + k - 1)*t + 2*pi*mod(0.6180339887*k^2, 1)); + end + x = x + nz * (2e-5*10^(p.Ln/20) / sqrt(mean(nz.^2))); + end + x = double(single(x)); + OUT = Tonality_Aures1985(x, p.fs, 0, 0, 0); + p.sum_x = sum(x); + p.sum_x2 = sum(x.^2); + p.K = OUT.InstantaneousTonality; + p.time = OUT.time; + p.Kmean = OUT.Kmean; + cases = il_append(cases, p); + fprintf('%-14s Kmean = %.4f t.u. (%d frames)\n', p.name, p.Kmean, numel(p.K)); +end +il_save('tonality_sqat.mat', cases, sqatRoot); +end + +function cases = il_append(cases, p) +if isempty(cases) + cases = p; +else + cases(end+1) = p; +end +end + +function il_save(name, cases, sqatRoot) %#ok +matlab_version = version; %#ok +sqat_commit = il_sqat_commit(sqatRoot); %#ok +out = fullfile(fileparts(mfilename('fullpath')), name); +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/tonality_sqat.mat b/tests/golden/tonality_sqat.mat new file mode 100644 index 0000000000000000000000000000000000000000..cbd2974536da48987c0124afba9be5176e89c288 GIT binary patch literal 1053 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X01nwjV*I2WZRmZYXAcz+3Ya0 z#@{3J?tE}4Th7e<@6VCeCGiZW*q>?!Gc#8lTav!-&moIh`s|_mQYtpK$Nbg)H0g`R z$?50A=4bBy^?R?a)n6HQ+Z&%%H(AR{e0Z)N^L6P-hO5)Lf4F~ss&Hz0-2eZd-YZA{ zKT@ab^(l4Qr-w0a8^gct$Vk(8yM}ckd-c9GfsQ*aH#Y96;mkQM#lN&{jf_}1-?l%d zET2_3Z2x+oE#_}I)BKm&Yv=fs)&1PBvD{|I-?_FMza*TWe&St8!h!kocTD`k?eKGf zO(@Tq8ErQw#OdjNRuyJ0`jyQ1>@=^=;;&C1oSe?vqg#FQYJI}H66NIN^i8k5pEo_S zjGQU|rgJIK2`*57w}brMermsK{s7hav-+L|&w3WDZ1h%o#jh_1O7yLZ zzJ~8;-y63+%d5;z_fGYpQ$dfq{%-!);$1s8{LTHpuR~3~?MaG%Gf!&r%z~r$o;dbi zyuU#E+oI=pdw0x_m;CZ`%fo~EUp;SzzhRWzBs>w6&V7IwlFmCAlFSsIF(iqxaU`^U zl5C8S@%qeIn3v?~mlRYJnH5-=Rwu{Y%n)+#gwJ`8^A6{4+&FaR!T~2o4wS0J6DVhEw5Ubo;-2n kj>W}G7w%k8GCX?Y#)(TNmPa)%DrxF!9%W*1O%gZ*09(M?xc~qF literal 0 HcmV?d00001 diff --git a/tests/test_psychoacoustics.py b/tests/test_psychoacoustics.py index a1bb3d3..b06ed26 100644 --- a/tests/test_psychoacoustics.py +++ b/tests/test_psychoacoustics.py @@ -5,8 +5,11 @@ standard, parity with SQAT (https://github.com/ggrecow/SQAT) on golden data, or the rejection of an invalid input. -The SQAT golden data come from golden/make_loudness_sqat.m and -golden/make_roughness_sqat.m. Those scripts store the parameters of short +The SQAT golden data come from golden/make_loudness_sqat.m, +golden/make_roughness_sqat.m and golden/make_tonality_sqat.m (SQAT's +tonality is GPL-3.0 and serves only as an oracle; pytta's code for it is +written from the papers). Those scripts +store the parameters of short synthetic signals; the helpers below rebuild the same samples (rounded to float32, as SignalObj stores them) and check them against the sums stored with the outputs. @@ -65,6 +68,22 @@ def roughness_signal(p): return x.astype(np.float32) +def tonality_signal(p): + """Rebuild a tonality case of make_tonality_sqat.m.""" + t = np.arange(int(round(p.dur * p.fs))) / p.fs + x = np.zeros(t.size) + for f in np.atleast_1d(p.f): + x = x + np.sqrt(2) * P0 * 10 ** (p.L / 20) * np.sin(2 * np.pi * f * t) + if p.kind == 'tone+noise': + K = int(p.f2 - p.f1 + 1) + nz = np.zeros(t.size) + for k in range(1, K + 1): + nz = nz + np.sin(2 * np.pi * (p.f1 + k - 1) * t + + 2 * np.pi * np.mod(0.6180339887 * k ** 2, 1)) + x = x + nz * (P0 * 10 ** (p.Ln / 20) / np.sqrt(np.mean(nz ** 2))) + return x.astype(np.float32) + + def same_samples(tc, p, x): """The rebuilt float32 signal is the one MATLAB used.""" x = x.astype(float) @@ -292,5 +311,61 @@ def test_rejects_invalid_input(self): pa.roughness_daniel1997(short) # shorter than one 200 ms window +class TestTonality(unittest.TestCase): + + def test_reference_tone_is_1_tu(self): + # Aures (1985), 3.1: a 1 kHz tone at 60 dB has a tonality of 1. + fs = 48000 + t = np.arange(fs) / fs + x = np.sqrt(2) * P0 * 1000 * np.sin(2 * np.pi * 1000 * t) + out = pa.tonality_aures1985(pytta.SignalObj(x, 'time', fs)) + print(f'\n Kmean = {out["Kmean"]:.6f} t.u., diff {out["Kmean"] - 1:+.1e}') + self.assertAlmostEqual(out['Kmean'], 1, delta=1e-5) + np.testing.assert_allclose(out['WGr'], 1, atol=1e-6) # no noise left + + def test_level_excess_weighting(self): + # A lone tone is masked by the threshold only, so its SPL excess is + # L - LTH(f) (Terhardt et al. 1982, Eqs. 4 and 8). With w2 and WGr + # unchanged, K(30 dB) / K(60 dB) = (w3'(30 - LTH) / w3'(60 - LTH))^0.29, + # w3' = 1 - exp(-LX / 15 dB) (Aures 1985, Eqs. 10, 11 and 13). + fs = 48000 + t = np.arange(fs) / fs + x = np.sqrt(2) * P0 * 10 ** (30 / 20) * np.sin(2 * np.pi * 1000 * t) + out = pa.tonality_aures1985(pytta.SignalObj(x, 'time', fs)) + lth = 3.64 - 6.5 * np.exp(-0.6 * (1 - 3.3) ** 2) + 1e-3 + w3 = 1 - np.exp(-(np.array([30, 60]) - lth) / 15) + expected = (w3[0] / w3[1]) ** 0.29 + print(f'\n K(30 dB) = {out["Kmean"]:.6f}, expected {expected:.6f}, ' + f'diff {out["Kmean"] - expected:+.1e}') + self.assertAlmostEqual(out['Kmean'], expected, delta=1e-5) + + def test_matches_sqat(self): + # Tones: 500 Hz and 4 kHz (w2), 30 dB (w3), two tones 1 and 1.2 kHz + # (mutual masking in the SPL excess) agree to 1e-4 t.u. A 1 kHz tone + # at 60 dB over a noise from 200 Hz to 5 kHz at 60 and 50 dB (masking + # by noise and WGr): within 0.015 t.u. (the noise loudness comes + # from the band levels of the spectrum here, from a resynthesised + # signal in SQAT). SQAT's frames differ, so the means are compared. + for p in golden_cases('tonality_sqat.mat'): + x = tonality_signal(p) + same_samples(self, p, x) + out = pa.tonality_aures1985(pytta.SignalObj(x, 'time', p.fs)) + d = out['Kmean'] - p.Kmean + print(f'\n {p.name}: Kmean = {out["Kmean"]:.5f}, SQAT {p.Kmean:.5f}, diff {d:+.1e}') + self.assertLess(abs(d), 0.015 if p.kind == 'tone+noise' else 1e-4) + + def test_rejects_invalid_input(self): + fs = 48000 + one = pytta.SignalObj(np.zeros(fs), 'time', fs) + with self.assertRaises(ValueError): + pa.tonality_aures1985(one, field='near') + two = pytta.SignalObj(np.zeros((fs, 2)), 'time', fs) + with self.assertRaises(ValueError): + pa.tonality_aures1985(two) + short = pytta.SignalObj(np.zeros(fs // 10), 'time', fs) + with self.assertRaises(ValueError): + pa.tonality_aures1985(short) # shorter than one 160 ms frame + + if __name__ == '__main__': unittest.main() From 40974867983d44535a460d02bac90013d3fce863 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 23:46:10 +0100 Subject: [PATCH 3/4] [FEATURE] add the fluctuation strength of Osses et al. to pytta.psychoacoustics fluctuation_strength_osses2016(signal, time_skip=0.0) gives FS(t) every 0.2 s over frames of 2 s, the specific fluctuation strength over 47 half-Bark channels and the mean FSmean, in vacil, after Osses, Garcia and Kohlrausch (2016): the roughness model of Daniel and Weber with frames of 2 s, a modulation weighting H(fmod), the generalised modulation depth compressed above 0.7 and FS = C_FS sum (m*)^1.7 |k k|^1.7 g(z). It is written from the paper. SQAT's FluctuationStrength_Osses2016 and its H(fmod) data are GPL-3.0 and served only as the test oracle. The paper does not give the filter type of H(fmod), so its cut-offs (1.75 and 13.5 Hz) were fitted to the AM tones of the paper's Table 1. The result differs from SQAT, as the docstring states: within 1 % for AM tones at 4 Hz, 4 % at 16 Hz, 14 % above for an FM tone, and about half of SQAT's value for AM at 1 Hz. Tests check the unit (1 vacil for the reference AM tone), the AM tones of Table 1 (rms deviation 0.09 vacil), the comparison with SQAT on golden data from tests/golden/make_fluctuation_sqat.m with the tolerances above, and the rejection of invalid input. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/index.rst | 4 +- pytta/__init__.py | 4 +- pytta/psychoacoustics.py | 168 ++++++++++++++++++++++++++- tests/golden/fluctuation_sqat.mat | Bin 0 -> 5407 bytes tests/golden/make_fluctuation_sqat.m | 94 +++++++++++++++ tests/test_psychoacoustics.py | 89 +++++++++++++- 6 files changed, 347 insertions(+), 12 deletions(-) create mode 100644 tests/golden/fluctuation_sqat.mat create mode 100644 tests/golden/make_fluctuation_sqat.m diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index acf770c..b86ef09 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -59,8 +59,8 @@ Modules: * pytta.psychoacoustics: loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel - and Weber 1997) and tonality (Aures 1985) of a sound pressure - signal; + and Weber 1997), fluctuation strength (Osses et al. 2016) and + tonality (Aures 1985) of a sound pressure signal; .. toctree:: classes/index diff --git a/pytta/__init__.py b/pytta/__init__.py index 6f4b03e..e7d3c16 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -74,8 +74,8 @@ * pytta.psychoacoustics: loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel - and Weber 1997) and tonality (Aures 1985) of a sound pressure - signal; + and Weber 1997), fluctuation strength (Osses et al. 2016) and + tonality (Aures 1985) of a sound pressure signal; You can find usage examples in the toolbox's examples folder. diff --git a/pytta/psychoacoustics.py b/pytta/psychoacoustics.py index 8f12817..363da69 100644 --- a/pytta/psychoacoustics.py +++ b/pytta/psychoacoustics.py @@ -9,12 +9,15 @@ >>> pytta.psychoacoustics.loudness_zwicker(signal) >>> pytta.psychoacoustics.sharpness_din45692(signal) >>> pytta.psychoacoustics.roughness_daniel1997(signal) + >>> pytta.psychoacoustics.fluctuation_strength_osses2016(signal) >>> pytta.psychoacoustics.tonality_aures1985(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); +* fluctuation_strength_osses2016: fluctuation strength after Osses, Garcia + and Kohlrausch (2016); * tonality_aures1985: tonality after Aures (1985), tonal components after Terhardt, Stoll and Seewann (1982). @@ -56,10 +59,12 @@ sharpness_din45692 is written from the text of DIN 45692:2009 (Formula (1)). -tonality_aures1985 is written from the papers cited below. SQAT's -Tonality_Aures1985 is GPL-3.0; it was run only as a test oracle, and of its -file only the help text (calling convention) was read. Where the papers -leave a choice open, the function docstring says what was chosen. +fluctuation_strength_osses2016 and tonality_aures1985 are written from the +papers cited below. SQAT's FluctuationStrength_Osses2016 and +Tonality_Aures1985 are GPL-3.0; they were run only as test oracles, and of +their files only the help text (calling convention) was read. The fluctuation strength reuses the Terhardt excitation +and the a0 factor of the roughness port (MIT). Where the papers leave a +choice open, the function docstrings say what was chosen. References ---------- @@ -74,6 +79,11 @@ Daniel, P., & Weber, R. (1997). Psychoacoustical roughness: implementation of an optimized model. Acustica, 83, 113-123. +Osses Vecchi, A., Garcia Leon, R., & Kohlrausch, A. (2016). Modelling the +sensation of fluctuation strength. Proceedings of Meetings on Acoustics, 28, +050005 (22nd International Congress on Acoustics, Buenos Aires). +doi:10.1121/2.0000410 + Aures, W. (1985). Berechnungsverfahren fuer den sensorischen Wohlklang beliebiger Schallsignale. Acustica, 59, 130-141. @@ -158,7 +168,7 @@ import scipy.signal as ss __all__ = ['loudness_zwicker', 'sharpness_din45692', 'roughness_daniel1997', - 'tonality_aures1985'] + 'fluctuation_strength_osses2016', 'tonality_aures1985'] _FIELDS = ('free', 'diffuse') _METHODS = ('stationary', 'time_varying') @@ -992,6 +1002,154 @@ def roughness_daniel1997(signal, time_skip=0.0): 'Rmean': float(np.mean(R[_skip_index(time, time_skip):]))} +# --------------------------------------------------------------------------- +# Fluctuation strength, Osses, Garcia and Kohlrausch (2016) +# --------------------------------------------------------------------------- + +_FS_WINDOW = 2.0 # s, analysis frame (Osses et al. 2016, 2.1) +_FS_HOP = 0.2 # s, 90 % overlap +_FS_RAMP = 0.05 # s, raised-cosine ramps at both ends of each frame +# H(fmod): cascade of a 2nd-order high-pass and a 4th-order low-pass (the +# orders of the paper, 3.1), Butterworth magnitudes. The paper gives the +# passband as 3.1 to 12 Hz; read as Butterworth cut-off frequencies those +# values give 0.35 vacil at fmod = 2 Hz and 1.6 at 8 Hz for the AM tones of +# its Table 1 (Fastl and Zwicker: 0.84 and 1.30). The cut-offs below are fitted +# to that table, as the paper fits its filter (least squares over the six +# fmod, C_FS re-set each time): rms deviation 0.09 vacil, and C_FS = 0.2480 +# against the paper's 0.2490. (Orders 1 and 2 would fit the table better, +# 0.03 vacil, with C_FS = 0.310.) +_FS_HP = 1.75 # Hz +_FS_LP = 13.5 # Hz +_FS_KNEE = 0.7 # compression of m* above this value, ratio 3:1 (2.1.3) +_FS_P = 1.7 # exponents pm = pk (3.1) +# C_FS of Eq. (1), set (as the paper does) so that the reference signal, a +# 1 kHz tone at 60 dB, 100 % amplitude modulated at 4 Hz, gives 1 vacil (mean +# over the 11 frames of 4 s of it). The paper reports 0.2490. +_C_FS = 0.248014 + + +def _fs_weights(N, fs): + """|H(fmod)| on the FFT bins of N samples and g(z) of the 47 channels.""" + f = np.abs(np.fft.fftfreq(N, 1 / fs)) + with np.errstate(divide='ignore'): + hp = np.where(f > 0, 1 / np.sqrt(1 + (_FS_HP / f) ** 4), 0.0) + lp = 1 / np.sqrt(1 + (f / _FS_LP) ** 8) + z = np.arange(1, _CHNO + 1) / 2 + g = np.interp(z, [15, 23.5], [1, 0.5]) + return hp * lp, g + + +def _specific_fluctuation(ei, H, g): + """Terms (m*)^pm |k_(i-2) k_i|^pk g(z_i) of Eq. (1), one per channel.""" + env = np.abs(ei) + h0 = np.mean(env, axis=1) + hbp = np.real(np.fft.ifft(np.fft.fft(env, axis=1) * H, axis=1)) + rms = np.sqrt(np.mean(hbp ** 2, axis=1)) + m = np.zeros(_CHNO) + pos = h0 > 0 + m[pos] = rms[pos] / h0[pos] + over = m > _FS_KNEE + m[over] = _FS_KNEE + (m[over] - _FS_KNEE) / 3 + hc = hbp - np.mean(hbp, axis=1, keepdims=True) + norm = np.sqrt(np.sum(hc ** 2, axis=1)) + k = np.zeros(_CHNO - 2) # k between channels i and i + 2 + for i in range(_CHNO - 2): + if norm[i] > 0 and norm[i + 2] > 0: + k[i] = np.dot(hc[i], hc[i + 2]) / (norm[i] * norm[i + 2]) + kk = np.ones(_CHNO) + kk[:-2] *= k # k_i: with the channel 1 Bark above + kk[2:] *= k # k_(i-2): with the channel 1 Bark below + return m ** _FS_P * np.abs(kk) ** _FS_P * g + + +def fluctuation_strength_osses2016(signal, time_skip=0.0): + """ + Fluctuation strength after Osses, Garcia and Kohlrausch (2016). + + Written from the paper; SQAT's FluctuationStrength_Osses2016 is GPL-3.0 + and served only as a test oracle. The model is the roughness model of + Daniel and Weber with these changes: frames of 2 s with 90 % overlap, + gated with 50 ms raised-cosine ramps; the a0 transmission factor and + the 47 half-Bark channels with Terhardt's slopes of + roughness_daniel1997; a modulation weighting H(fmod), a 2nd-order + high-pass times a 4th-order low-pass applied to the spectrum of each + channel envelope; the generalised modulation depth m* (rms of the + weighted envelope over its mean) compressed 3:1 above 0.7; and + + FS = C_FS * sum_i (m*_i)^1.7 |k_(i-2) k_i|^1.7 g(z_i), + + with k the normalised cross covariance of the weighted envelopes of + channels 1 Bark apart. 1 vacil: 1 kHz tone at 60 dB, 100 % amplitude + modulated at 4 Hz. Signals at rates other than 40960, 44100 and 48000 + Hz are resampled to 48 kHz first. + + The paper leaves open the filter type of H(fmod), the shape of g(z) + and its a0 curve. Here: Butterworth magnitudes with cut-offs 1.75 and + 13.5 Hz, fitted to the AM tones of the paper's Table 1 as the paper + fits its own filter (the 3.1 and 12 Hz it quotes, taken as Butterworth + cut-offs, miss that table by up to 0.5 vacil); g(z) = 1 up to 15 Bark, + then linear to 0.5 at 23.5 Bark ("reducing g(z_i) gradually"); a0 for + the free field after Fastl and Zwicker (2007), Fig. 8.18 (SQAT's option + a0_type = 'fastl2007'; its default a0 is flat up to 5 kHz). Against + SQAT: within 1 % for AM tones at 4 Hz, 4 % at 16 Hz, 14 % above for an + FM tone, and about half of SQAT's value at fmod = 1 Hz. + + Parameters (default), (type): + ----------------------------- + + * signal (), (SignalObj): + One channel, sound pressure in Pa, at least 2 s long; + + * time_skip (0.0), (float): + Seconds excluded from FSmean at the start. + + Returns: + -------- + + dict, nt frames every 0.2 s: + + * 'FS': fluctuation strength FS(t), vacil, shape (nt,); + * 'specific': specific fluctuation strength, vacil/Bark, + shape (nt, 47) (FS = 0.5 * sum over the channels); + * 'time': start of each frame, s, shape (nt,); + * 'bark': critical band rate of 'specific', 0.5 to 23.5 Bark, + shape (47,); + * 'FSmean': mean of FS(t) from the frame closest to time_skip + on, vacil (float). + + """ + x, fs = _mono(signal) + if fs not in (44100, 40960, 48000): + x = _resample(x, 48000, fs) + fs = 48000 + N = _mround(fs * _FS_WINDOW) + hop = _mround(fs * _FS_HOP) + if x.size < N: + raise ValueError(f'the signal ({x.size / fs:.3f} s) is shorter than ' + f'one frame ({_FS_WINDOW} s).') + nwin = (x.size - N) // hop + 1 + prm = _excitation_params(N, fs) + H, g = _fs_weights(N, fs) + nr = _mround(fs * _FS_RAMP) + window = np.ones(N) + ramp = 0.5 - 0.5 * np.cos(np.pi * np.arange(nr) / nr) + window[:nr] = ramp + window[N - nr:] = ramp[::-1] + 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)) + spec = np.zeros((nwin, _CHNO)) + time = np.arange(nwin) * hop / fs + for w in range(nwin): + start = w * hop + frame = np.fft.fft(x[start:start + N] * window * amp_cal) * prm['a0'] + ei, _ = _terhardt_excitation(frame, prm) + spec[w] = _C_FS / _R_DZ * _specific_fluctuation(ei, H, g) + F = _R_DZ * np.sum(spec, axis=1) + return {'FS': F, 'specific': spec, 'time': time, + 'bark': np.arange(1, _CHNO + 1) / 2, + 'FSmean': float(np.mean(F[_skip_index(time, time_skip):]))} + + # --------------------------------------------------------------------------- # Tonality, Aures (1985), tonal components after Terhardt et al. (1982) # --------------------------------------------------------------------------- diff --git a/tests/golden/fluctuation_sqat.mat b/tests/golden/fluctuation_sqat.mat new file mode 100644 index 0000000000000000000000000000000000000000..964c27865a5c6dc627d7a142230239fbc0abeafe GIT binary patch literal 5407 zcma)=Ra6v=x5ep3O1eR5M_TFbl9C#vqD3)uZ7UtoYDpQQ_!&&&068UW2lY{&>@8!d$9$RIl$AxC#LG?lQ^fB&IU9Lk~lFX~pF zx-7SAu%;z4GOsKGm!H~PSi|al*?mSenyI0r{YVm?%dEUcBPuS0^Q(E%Q^PDH>D->y zg7Xi%;^obKD+csfh~F!Ynn4bN%|)}1QMWWPw&M<3I+3?@JKlFiHI(hw{ekG0YXt=n z;=&KVhZ^B8xlFE2Gv3;ox+W_4Br4QX?r6`ro8gY50Um?=!fJ7jw!h*lrdmq9;<}wL z%AN#PKGMOJ((MQj7KF1j{T#J-%N*q|-3uOV)%Kuxt~~#IK9`H=qC{N0Z;y6ZynKJ$ z!1VFvk|?+%3A02~X_W7v11)Ay+*StjzD1YVt-*xr+qgcilZ`R96-OOijyK8hDyvbq++s~ogMmNgRsNME)}`0XF%%(XY335*W`!4MU4#M+yA!PV>oLM!<{_SZ!4CgX|%* zdYIQFtX7do|1hFJt;4$85GIX_AD%6TL2RCYYJDVwb^Wy7;2@GS4f3 zpy-6S$u+dlW4hT6HbOrC;-_$>C?MmE(Ua3?WpQS=gm}25cU_Yh7YOGa8 zGLHnX*Z0BPNtlM+=s9^=C*p11N0&ujPDyFp?D>&L%j%cSPk3HrY_9Af)2A2ELAdSw$VoPveu5lyQrKvm4;R3o+@1VMavgACncLY9&BiF z;J&gsdGj6{BO0@+kIXhOANBrrC>#*R->d2XcJeUbbeL&Ml2Gd_L%aD_{F=4VO8FR? zcivFMRN+Xy-rTR9yJtNDRAS${eEYthMtixF7PTnO3b8a9xtVJ3lJIS+OFhbpZ8a@_ zJ_6-(KpXw_p)lFwW+DM5koQIb2RA2w&s%`2!Y~0Z-%rkHqx&gP)(+s;A1%ME`~b3X zDs?PczM-}@kj>mC1#Vd1m^ZyE;s3;+I$~Z+Vn7R~%c-niNpG@PFMrCsY`bO1J!|mo zZBb#eRYC)wK4PefQ8B=zVE3p#G1b~%AtYP;{m3n^jB-6-I6Dk{NQ=P>$-ug!Y{2Yh z2kCta-K9@&SFZ-=t+$3|kNh|Lsdk5&g3rQk(}{frp5y`#`zU@IUx218(JmjD9j!PF63R0-$hKd|TYHNItoOcA zvfgD)Rmwyb@4VPDm4FLOf7^c=47j;X@2kl(MZQ5*$SJccJoWcKk*k_vFLm;F3W{Gc zWK&VG5UU;il0%%8ja;t2$?QD&GWz<@#hDte>(p>N8WiW{+$(q9 zsDQ_y(+!W`SjOaB_suYfPM&|tfb4~bL-41g#@Ad3Pfw%9ld3X$>Bn{|s>du6#{Dxf z)7`Yf#{N;@`-(RN;H?nt{|3v}Tib{1Z89^R z?@|KC|1uwcXsh^>hkQMuF(%otdUaLl=$BsFjPL*G&+zGOWx=sV4O%kNYy_YFeM~4u?@ZD9OmGySv2Zm}j6dvQUxE!dnRyg`q+H0_sguzc{8MRlm z%D=h`I>j?zO=^Cf)?vLOKz-B(qZM5t+mEzMmeo4P2Osc8$o+uS7JI(pY&KHLe?}3# zSZQ?dbNZD7K76Sy7p1-G-*I8C!a@1$+jVC4`xVt-Cj+I@&@sJ0}ynF2v z+L*-CfL+@%6r$z-?g{tFLw6V~jRRKF7pN*}K@i%au8hHMA!lu?R}MVp;N;UthWx?U^FH|lKKP1S@w(_TAqV%nz1 zQlPdfIIXy%HhQ^Q3u>Rx)DWqdF3e05V^?;sAe#KA!SVdCnTwIVh&cDTY<~TsHb9=R zr+p88^jdFPkF{s$tUN#Y2J>F9+VrxmDrvIupL)Uxa1f{T=F2dI%ilAY2L}#(jd$2- zSPLA3c}jsQLa={5;XhIY!gXuB@9IBWu>l)9+}eZ0fRpe#z_~%nVF2h+NePu<@9sPM z;o1%VEa=|35OQ}^> zXg9M@)I%EE@p-Yy{ELvi>dD`-$(-pb_tt|h^3sc(UxCWTzk9FjlQpTqE1!C#8xAc< zLGQBDg3w)N#Nq$wjRh}T#xjQYD@q_KDFxE`5-9}WjBQj(9j#h4O2tz!PeH|SlhW|- zx~|w;2Aj5)iW@UY2Q!$>Z29BBi@1LW}n?jeqGYRYW-b|OQE4RmL zI`Fox|HDdM!N-LkQD~k0F_wjyxaq53!^tu{g6groO<>nCiVuNu)Z)QB37ZMq^=J26 zLtYqQRZ>`=8R`YJs`5rj9B!-mILb`Knc#dgnT{v5?2KU070wSv;#8D=-~s7{im%kq>B zZ$NA$y!HitBb8^RG08vRE<6J%jgz8SB9HgofqCmA`4gYfiNjH`q$A=$=W$l6r*^?E zUP;PWv1Yi}7;IakQ7{trL4al_QX`uj zY0>9PCUma**0J$hN1x6&dW$q(FQpXE9|P>HW-Q4rzH$Fn>Fw)kI>qQORX6l9$E{t)?3mo z$27;QB8Lde>v)jdt*C-eR<_ip6$u^3VGD`q4?ghvP5TJ?mGD0^@^XC9uA%xJvYam> z1Z(I=spi(;WTg3ir&XB}^zjN4#+7Xv%4@|z9iKQ7Lr3Ijmn9Q;UI8gZL5+bjOsQcz z9Ut@F&|e}C^lIgWfQrsv5CU~nH1%ccbqwMaG1Q>C=T}i0};x59zo_ue%d#Y zM&suGYz}3hSHaDywBdKh$5sE{^XXUG!8=xtym8v5i zYfD%>?s*@_8)j>!(%5rds9!zbcEr68S>RE+%M&Hbe>-Ci0J$;NcJ0pG&Z!GqSk~nC z8Ecr%Q}uAkgT5~FFAH9pAs@-lq}j1d7a-)|$}94w{P$H9)kiMb&NzS03$F0AwJ3q)$P8Mu;JT0r-VjB;@ys~B*_YlIVgkD#I2{_r2d~}KH^)^vJUe zdj|bfx5~EMbvvdp!%4NB z$z@L(;IedVEK$EfIkU4TtW^{gqKnDX_I;=y#^2@n7{t<*M=ho^9mCe+#te7=fzw;n zy`N-KaPj-p%xI`-b6lu%+I5BgACLYrKq4?e)IqqO0_2((f`TATf0`f_Tg++@AFZl#_v55&~#(rNa92)$$$lu)^O&ypOc(RcNJUIv@{j;3THd3CFiahJ9ZoHr~b6C zH(}!#3Z~F7 zub2)3Ntcpx#iF5|KynAdGOpLwbgf?<_D|n(I#@yTxfIr>5_Z;PXIJL!)-J&P)bZI>^>X7!d*)SFEh5lWKI5&!(b zZre7GMZoAi?tUVetilIOa5IU0H)$uj?tyy9sZ$+w;ZT%*MrUPKav%uhQDUL5f5QRu zC+}5a1e;4?=1i;H;=vi9*xPJ~JCo0EaRz(+hd$j`McBPz3pv~6MxJuu^Caotb48+) zQlEqn!Gg0z^RulXh{o7VtR;^kIeF8|6(G4;K>OTHNyqG}qQwJ_Iii!ds*A04DGdSD z>hbF2laGx=yKK^$gh?E0Al$pXREP=9)?9tVCEbgCK3R)cJa6m0N`--`934fG?_b8h zM6R~Tgs8-1gsM+!+62*}Qovx?Fviez_SP-|3sBnm{&HY*J;4^6&)bo=ukfstxsmj) zm{rAz#AkRLc3fQF`?X>M29|Q1c7zm8ZdSEJbo7?0+*)Vzr&pF41^S~dtEpZ} zqcUzkQZYKx^_7~`50yo`lt-3Rk;SXBTkTgT;+}!cRuPn3vnM#rC=u9WT(#iH&r@j? zSMKJ^7C&OkHa-_I0O<*Say#<{QgnhZd6aV;beHo!o>1=IIop3fea8Pz1a#ZY**TrS zRi1B;z?FrJ!a3!XJ;9dC@Qg4xz}3Li)gs2u$2{00j)ypjz-kR{JflCOKC`m&&-d9O zO$|Xfd +matlab_version = version; %#ok +sqat_commit = il_sqat_commit(sqatRoot); %#ok +out = fullfile(fileparts(mfilename('fullpath')), name); +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/test_psychoacoustics.py b/tests/test_psychoacoustics.py index b06ed26..6498e8c 100644 --- a/tests/test_psychoacoustics.py +++ b/tests/test_psychoacoustics.py @@ -6,9 +6,10 @@ or the rejection of an invalid input. The SQAT golden data come from golden/make_loudness_sqat.m, -golden/make_roughness_sqat.m and golden/make_tonality_sqat.m (SQAT's -tonality is GPL-3.0 and serves only as an oracle; pytta's code for it is -written from the papers). Those scripts +golden/make_roughness_sqat.m, golden/make_fluctuation_sqat.m and +golden/make_tonality_sqat.m (SQAT's fluctuation strength and tonality are +GPL-3.0 and serve only as an oracle; pytta's code for them is written from +the papers). Those scripts store the parameters of short synthetic signals; the helpers below rebuild the same samples (rounded to float32, as SignalObj stores them) and check them against the sums stored @@ -68,6 +69,18 @@ def roughness_signal(p): return x.astype(np.float32) +def fluctuation_signal(p): + """Rebuild a fluctuation case of make_fluctuation_sqat.m.""" + t = np.arange(int(round(p.dur * p.fs))) / p.fs + if p.kind == 'am': + A = np.sqrt(2) * P0 * 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) + else: + x = (np.sqrt(2) * P0 * 10 ** (p.L / 20) + * np.sin(2 * np.pi * p.fc * t - (p.df / p.fm) * np.cos(2 * np.pi * p.fm * t))) + return x.astype(np.float32) + + def tonality_signal(p): """Rebuild a tonality case of make_tonality_sqat.m.""" t = np.arange(int(round(p.dur * p.fs))) / p.fs @@ -311,6 +324,76 @@ def test_rejects_invalid_input(self): pa.roughness_daniel1997(short) # shorter than one 200 ms window +class TestFluctuationStrength(unittest.TestCase): + + def test_reference_am_tone_is_1_vacil(self): + # Definition of the vacil: 1 kHz tone, 100 % amplitude modulated at + # 4 Hz, 60 dB (rms of the modulated signal). C_FS is set on 4 s of + # this signal, so the mean of its 11 frames is 1 to the 6 digits of + # C_FS. + fs = 48000 + t = np.arange(4 * fs) / fs + A = np.sqrt(2) * P0 * 1000 / np.sqrt(1.5) + x = A * (1 + np.sin(2 * np.pi * 4 * t)) * np.sin(2 * np.pi * 1000 * t) + out = pa.fluctuation_strength_osses2016(pytta.SignalObj(x, 'time', fs)) + print(f'\n FSmean = {out["FSmean"]:.6f} vacil, diff {out["FSmean"] - 1:+.1e}') + self.assertAlmostEqual(out['FSmean'], 1, delta=1e-5) + self.assertEqual(out['specific'].shape, (11, 47)) + np.testing.assert_allclose(out['FS'], 0.5 * np.sum(out['specific'], axis=1)) + + def test_am_tones_follow_fastl_and_zwicker(self): + # Osses et al. (2016), Table 1: FS of 1 kHz AM tones (m = 1) at 70 dB + # after Fastl and Zwicker, the data the paper fits H(fmod) to and + # that pytta's cut-offs are fitted to (rms 0.089 vacil, worst 1 Hz, + # -0.18 vacil). Pins the band-pass shape: 0.1 vacil rms, 0.2 worst. + fs = 48000 + t = np.arange(4 * fs) / fs + A = np.sqrt(2) * P0 * 10 ** (70 / 20) / np.sqrt(1.5) + table = {1: 0.39, 2: 0.84, 4: 1.25, 8: 1.30, 16: 0.36, 32: 0.06} + dev = [] + for fm, ref in table.items(): + x = A * (1 + np.sin(2 * np.pi * fm * t)) * np.sin(2 * np.pi * 1000 * t) + F = pa.fluctuation_strength_osses2016(pytta.SignalObj(x, 'time', fs))['FSmean'] + dev.append(F - ref) + print(f'\n fmod {fm:2d} Hz: FS = {F:.3f} vacil, table {ref:.2f}, diff {F - ref:+.3f}', + end='') + rms = np.sqrt(np.mean(np.square(dev))) + print(f'\n rms deviation {rms:.3f} vacil') + self.assertLess(rms, 0.1) + self.assertLess(np.max(np.abs(dev)), 0.2) + + def test_matches_sqat(self): + # SQAT's H(fmod) is not described in the paper (it is loaded from + # files), so the two codes agree where H is flat: within 1 % at 4 Hz + # (48 and 44.1 kHz), 5 % at 16 Hz. The FM tone (1500 Hz, +-700 Hz, + # 4 Hz, 70 dB; Table 1 gives 2.00 vacil) is 14 % above SQAT; the + # paper notes that its model overestimates FM tones. At 1 Hz the two + # differ by 46 % (see the Fastl and Zwicker test), so it is left out. + # The frames start at the same samples (SQAT's time is 1 sample + # later). + tolerance = {'am1k_4hz': 0.01, 'am1k_4hz_441': 0.01, 'am1k_16hz': 0.05, + 'fm1k5_4hz': 0.15} + for p in golden_cases('fluctuation_sqat.mat'): + x = fluctuation_signal(p) + same_samples(self, p, x) + out = pa.fluctuation_strength_osses2016(pytta.SignalObj(x, 'time', p.fs)) + self.assertEqual(out['FS'].shape, p.FS.shape) + np.testing.assert_allclose(out['time'], p.time - 1 / p.fs, atol=1e-9) + rel = out['FSmean'] / p.FSmean - 1 + print(f'\n {p.name}: FSmean = {out["FSmean"]:.4f}, SQAT {p.FSmean:.4f}, ' + f'rel diff {rel:+.4f}') + self.assertLess(abs(rel), tolerance[p.name]) + + def test_rejects_invalid_input(self): + fs = 48000 + two = pytta.SignalObj(np.zeros((3 * fs, 2)), 'time', fs) + with self.assertRaises(ValueError): + pa.fluctuation_strength_osses2016(two) + short = pytta.SignalObj(np.zeros(fs), 'time', fs) + with self.assertRaises(ValueError): + pa.fluctuation_strength_osses2016(short) # shorter than one 2 s frame + + class TestTonality(unittest.TestCase): def test_reference_tone_is_1_tu(self): From a2544059227da7f790ad36544152ec0eee756f2e Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 23:57:29 +0100 Subject: [PATCH 4/4] [FEATURE] add sharpness weightings, loudness level and critical band to pytta.psychoacoustics sharpness_din45692 gets weighting='din', 'aures' or 'bismarck': the weightings of Aures (which depends on the total loudness) and of von Bismarck, Formulae (B.1) and (B.2) of DIN 45692, Annex B, with the same k = 0.11, as SQAT and the MATLAB acousticSharpness that the ITA-Toolbox calls offer. loudness_level and loudness_from_level convert between sone and phon by ISO 532-1, Formulae (1) to (3), as ita_sone2phon. critical_band gives the critical band rate and bandwidth of Zwicker and Terhardt (1980), with the band edges, as ita_freq2cb (which squares the arctangent of the second term of the rate; the docstring says so). Tests check that von Bismarck equals DIN when no loudness lies above 15 Bark, parity of the three weightings with SQAT on golden data from tests/golden/make_sharpness_weights_sqat.m (within 3.6e-15 acum), the loudness level formulae and their inverse, and the critical band at 1 kHz and against the table of Zwicker. Co-Authored-By: Claude Fable 5.1 --- pytta/psychoacoustics.py | 96 ++++++++++++++++++--- tests/golden/make_sharpness_weights_sqat.m | 41 +++++++++ tests/golden/sharpness_weights_sqat.mat | Bin 0 -> 573 bytes tests/test_psychoacoustics.py | 80 ++++++++++++++++- 4 files changed, 204 insertions(+), 13 deletions(-) create mode 100644 tests/golden/make_sharpness_weights_sqat.m create mode 100644 tests/golden/sharpness_weights_sqat.mat diff --git a/pytta/psychoacoustics.py b/pytta/psychoacoustics.py index 363da69..9404008 100644 --- a/pytta/psychoacoustics.py +++ b/pytta/psychoacoustics.py @@ -11,6 +11,8 @@ >>> pytta.psychoacoustics.roughness_daniel1997(signal) >>> pytta.psychoacoustics.fluctuation_strength_osses2016(signal) >>> pytta.psychoacoustics.tonality_aures1985(signal) + >>> pytta.psychoacoustics.loudness_level(N) + >>> pytta.psychoacoustics.critical_band(f) * loudness_zwicker: Zwicker loudness, ISO 532-1:2017, stationary and time-varying sounds, free or diffuse field; @@ -87,6 +89,10 @@ Aures, W. (1985). Berechnungsverfahren fuer den sensorischen Wohlklang beliebiger Schallsignale. Acustica, 59, 130-141. +Zwicker, E., & Terhardt, E. (1980). Analytical expressions for critical-band +rate and critical bandwidth as a function of frequency. The Journal of the +Acoustical Society of America, 68(5), 1523-1525. doi:10.1121/1.385079 + Terhardt, E., Stoll, G., & Seewann, M. (1982). Algorithm for extraction of pitch and pitch salience from complex tonal signals. The Journal of the Acoustical Society of America, 71(3), 679-688. doi:10.1121/1.387544 @@ -168,7 +174,8 @@ import scipy.signal as ss __all__ = ['loudness_zwicker', 'sharpness_din45692', 'roughness_daniel1997', - 'fluctuation_strength_osses2016', 'tonality_aures1985'] + 'fluctuation_strength_osses2016', 'tonality_aures1985', + 'loudness_level', 'loudness_from_level', 'critical_band'] _FIELDS = ('free', 'diffuse') _METHODS = ('stationary', 'time_varying') @@ -547,6 +554,51 @@ def _loudness_level(N): return LN +def loudness_level(N): + """ + Loudness level LN (phon) of a loudness N (sone), ISO 532-1:2017, + Formulae (2) and (3): 40 + 33.22 lg N from 1 sone up, 40 (N + + 0.0005)^0.35 below (as ita_sone2phon). loudness_zwicker reports LN as + the program of Annex A does, 40 N^0.35 and at least 3 phon below 1 + sone, which differs from Formula (3) by 0.01 phon at 0.5 sone. + """ + N = np.asarray(N, dtype=float) + with np.errstate(divide='ignore'): + return np.where(N >= 1, 40 + 10 * np.log2(np.maximum(N, 1e-300)), + 40 * (np.maximum(N, 0) + 0.0005) ** 0.35) + + +def loudness_from_level(LN): + """ + Loudness N (sone) of a loudness level LN (phon), the inverse of + loudness_level: 2^((LN - 40) / 10) from 40 phon up (ISO 532-1, + Formula (1)), (LN / 40)^(1 / 0.35) - 0.0005 below. + """ + LN = np.asarray(LN, dtype=float) + return np.where(LN >= 40, 2 ** ((LN - 40) / 10), + (np.maximum(LN, 0) / 40) ** (1 / 0.35) - 0.0005) + + +def critical_band(f): + """ + Critical band rate and bandwidth at f (Hz), Zwicker and Terhardt (1980), + as ita_freq2cb: z = 13 arctan(0.76 f/kHz) + 3.5 arctan((f/7.5 kHz)^2) + Bark and df = 25 + 75 (1 + 1.4 (f/kHz)^2)^0.69 Hz, with the band edges + f_low f_high = f^2 and f_high - f_low = df. (ita_freq2cb squares the + arctangent of the second term, which gives 0.6 Bark less at 8 kHz.) + + Returns + ------- + dict + 'bark', 'bandwidth', 'f_low' and 'f_high', arrays like f. + """ + f = np.asarray(f, dtype=float) + df = 25 + 75 * (1 + 1.4 * (f / 1000) ** 2) ** 0.69 + f_low = 0.5 * (-df + np.sqrt(df ** 2 + 4 * f ** 2)) + return {'bark': _bark_terhardt(f), 'bandwidth': df, 'f_low': f_low, + 'f_high': f_low + df} + + def _temporal_weighting(N): """Step 10: two first-order low-passes (3.5 ms and 70 ms) combined.""" dec = _FS // _SR_LEVEL @@ -677,29 +729,46 @@ def loudness_zwicker(signal, field='free', method='stationary', time_skip=0.0): _K_SHARPNESS = 0.11 -def _sharpness(specific): +_WEIGHTINGS = ('din', 'aures', 'bismarck') + + +def _sharpness(specific, weighting='din'): """ - 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. + DIN 45692 Formula (1), or Formulae (B.1) and (B.2) of its Annex B, 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) + if weighting == 'din': + g = np.where(z <= 15.8, 1.0, 0.15 * np.exp(0.42 * (z - 15.8)) + 0.85) + elif weighting == 'bismarck': + g = np.where(z <= 15, 1.0, 0.2 * np.exp(0.308 * (z - 15)) + 0.8) + else: + # Aures: g_A(z) = 0.078 exp(0.171 z) / z * N / ln(0.05 N + 1), with + # N the total loudness (N / ln(0.05 N + 1) -> 20 as N -> 0) + N = np.atleast_1d(den * 0.1)[..., None] + with np.errstate(invalid='ignore', divide='ignore'): + ratio = np.where(N > 0, N / np.log(0.05 * N + 1), 20.0) + g = (0.078 * np.exp(0.171 * z) / z * ratio).reshape(np.shape(specific)) + num = np.sum(specific * g * z, 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): +def sharpness_din45692(signal, field='free', method='stationary', time_skip=0.0, + weighting='din'): """ 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. + should be reported as such. The weightings of Aures (S_A, which depends + on the total loudness N) and of von Bismarck (S_B), Formulae (B.1) and + (B.2) of its informative Annex B, use the same k; the standard asks to + label each result with its method. Parameters (default), (type): ----------------------------- @@ -715,7 +784,10 @@ def sharpness_din45692(signal, field='free', method='stationary', time_skip=0.0) * time_skip (0.0), (float): Passed to loudness_zwicker; for the time-varying method it also - starts the frames of S50. + starts the frames of S50; + + * weighting ('din'), (str): + 'din' (Formula (1)), 'aures' (B.1) or 'bismarck' (B.2). Returns: -------- @@ -733,8 +805,10 @@ def sharpness_din45692(signal, field='free', method='stationary', time_skip=0.0) single value DIN 45692 (clause 6) recommends for checks. """ + if weighting not in _WEIGHTINGS: + raise ValueError(f"weighting must be one of {_WEIGHTINGS}.") out = loudness_zwicker(signal, field=field, method=method, time_skip=time_skip) - S = _sharpness(out['specific']) + S = _sharpness(out['specific'], weighting) if method == 'stationary': return {'S': float(S)} tail = S[_skip_index(out['time'], time_skip):] diff --git a/tests/golden/make_sharpness_weights_sqat.m b/tests/golden/make_sharpness_weights_sqat.m new file mode 100644 index 0000000..1a566fe --- /dev/null +++ b/tests/golden/make_sharpness_weights_sqat.m @@ -0,0 +1,41 @@ +function make_sharpness_weights_sqat(sqatRoot) +% MAKE_SHARPNESS_WEIGHTS_SQAT Golden data for tests/test_psychoacoustics.py +% (sharpness with the weightings of Aures and von Bismarck, DIN 45692 Annex B). +% +% make_sharpness_weights_sqat(sqatRoot) runs Sharpness_DIN45692 of SQAT +% (https://github.com/ggrecow/SQAT) with the weightings 'DIN45692', +% 'aures' and 'bismarck', stationary free-field loudness, on band noises +% and writes sharpness_weights_sqat.mat next to this file. SQAT is only the +% test oracle (that file is GPL-3.0; pytta follows the text of the +% standard). +% +% Signals: multisine noise with components every 1 Hz from f1 to f2, +% phases 2*pi*mod(g*k^2, 1), g = 0.6180339887, k = 1..K, scaled to the +% rms level L (dB re 20 uPa), 1 s at 48 kHz, rounded to single precision. +% +% Usage: matlab -batch "make_sharpness_weights_sqat('/path/to/SQAT')" + +addpath(sqatRoot); +startup_SQAT([strip(sqatRoot, "right", filesep) filesep]); +defs = [920 1080 40; 920 1080 60; 920 1080 80; 3700 4400 60; 250 10000 60; 200 300 70]; +fs = 48000; t = (0:fs-1)' / fs; +cases = struct([]); +for c = 1:size(defs, 1) + p.f1 = defs(c,1); p.f2 = defs(c,2); p.L = defs(c,3); + K = p.f2 - p.f1 + 1; + x = zeros(size(t)); + for k = 1:K + x = x + sin(2*pi*(p.f1 + k - 1)*t + 2*pi*mod(0.6180339887*k^2, 1)); + end + x = double(single(x * (2e-5*10^(p.L/20) / sqrt(mean(x.^2))))); + p.sum_x = sum(x); + for w = {'DIN45692', 'aures', 'bismarck'} + OUT = Sharpness_DIN45692(x, fs, w{1}, 0, 1, 0, 0, 0); + p.(['S_' w{1}]) = OUT.Sharpness; + end + if isempty(cases), cases = p; else, cases(end+1) = p; end %#ok + fprintf('%5d-%5d Hz %2d dB: S %.4f, SA %.4f, SB %.4f\n', p.f1, p.f2, p.L, ... + p.S_DIN45692, p.S_aures, p.S_bismarck); +end +save(fullfile(fileparts(mfilename('fullpath')), 'sharpness_weights_sqat.mat'), 'cases', 'fs'); +end diff --git a/tests/golden/sharpness_weights_sqat.mat b/tests/golden/sharpness_weights_sqat.mat new file mode 100644 index 0000000000000000000000000000000000000000..75d7f33958c52510ee09b3624e737e0455c81593 GIT binary patch literal 573 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X01nwjV*I2WZRmZYXAcEk6rMF) zQaHf%)cBC0uyd3lGc)&uyc=f^Ffb^jG1Y^Vw=u$$XFNzsNMODs<8bB(1LK+iJ7s2< zK0&Q%b0X$M@uY1wC~!Vn+h}oi0wdd{o(Dei7TpVXn3(TlI>ggg{_%^kq3LE5CL>|* z?gMO0Gd*V>-P!4S@}`2EynwWU!kwdSM`I0mm~K6E1et1pE<0n+t=)^ zGU;vWiEIB=!{7898}DYeeE6-qaZiO%fc%W2{G?Ffn@{J}OseAd*XF+YQ0C*iyR2gJ zb(_8}*Pg#S-v&yxtr>W?#nUQmsZXJ0Fz= 1, 40 + 10 * np.log2(N), 40 * (N + 0.0005) ** 0.35) + back = pa.loudness_from_level(LN) + print(f'\n LN {np.round(LN, 3)}; round trip {np.abs(back - N).max():.1e} sone') + np.testing.assert_allclose(LN, expected, rtol=0, atol=1e-12) + self.assertEqual(float(pa.loudness_level(1)), 40) + self.assertAlmostEqual(float(pa.loudness_level(2)), 50, places=12) + np.testing.assert_allclose(back, N, rtol=1e-12) + + def test_critical_band_follows_zwicker_and_terhardt(self): + # At 1 kHz the bandwidth is 162.2 Hz with edges whose product is + # f^2 (ita_freq2cb gives 922.18 and 1084.39 Hz). The analytic rate + # stays within 0.2 Bark of the critical band table of Zwicker (lower + # edges of bands 1 to 24, at 1 to 24 Bark); the worst is 0.197 Bark + # at 4.4 kHz. + cb = pa.critical_band(1000) + print(f"\n 1 kHz: {cb['f_low']:.2f} to {cb['f_high']:.2f} Hz, " + f"{cb['bandwidth']:.4f} Hz, {cb['bark']:.4f} Bark", end='') + self.assertAlmostEqual(float(cb['bandwidth']), 162.2167, places=4) + self.assertAlmostEqual(float(cb['f_low'] * cb['f_high']), 1e6, places=6) + self.assertAlmostEqual(float(cb['f_low']), 922.1755, places=4) + edges = pa._BARK_RAW[1:, 1] + err = np.abs(pa.critical_band(edges)['bark'] - pa._BARK_RAW[1:, 0]).max() + print(f'; worst against the table {err:.2f} Bark') + self.assertLess(err, 0.2) + + class TestSharpness(unittest.TestCase): def test_din_reference_noise_is_1_acum(self): @@ -264,6 +309,37 @@ def test_time_varying_s50_of_reference_noise_is_1_acum(self): self.assertEqual(out['S'].shape, out['time'].shape) self.assertLessEqual(abs(out['S50'] - 1), 0.05) + def test_bismarck_equals_din_below_15_bark(self): + # g_B and g are both 1 up to 15 Bark (DIN 45692, (1) and (B.2)), so + # a sound with no specific loudness above 15 Bark has S_B = S. + x = multisine(920, 1080, 40) + sig = pytta.SignalObj(x, 'time', 48000) + spec = pa.loudness_zwicker(sig)['specific'] + S = pa.sharpness_din45692(sig)['S'] + SB = pa.sharpness_din45692(sig, weighting='bismarck')['S'] + print(f'\n N\' above 15 Bark: {spec[150:].max():.1e}; S {S:.6f}, S_B {SB:.6f}') + self.assertEqual(spec[150:].max(), 0) + self.assertAlmostEqual(S, SB, places=12) + + def test_weightings_match_sqat(self): + # Same float32 samples; DIN 45692, Aures (B.1) and von Bismarck (B.2) + # weightings on band noises from 40 to 80 dB. + d = sio.loadmat(os.path.join(GOLDEN, 'sharpness_weights_sqat.mat'), + squeeze_me=True, struct_as_record=False) + worst = 0 + for p in np.atleast_1d(d['cases']): + x = multisine(p.f1, p.f2, p.L) + self.assertAlmostEqual(np.sum(x, dtype=float), p.sum_x, delta=1e-9) + sig = pytta.SignalObj(x, 'time', 48000) + for w, ref in (('din', p.S_DIN45692), ('aures', p.S_aures), + ('bismarck', p.S_bismarck)): + S = pa.sharpness_din45692(sig, weighting=w)['S'] + worst = max(worst, abs(S - ref)) + print(f'\n {p.f1}-{p.f2} Hz {p.L} dB: S_A {pa.sharpness_din45692(sig, weighting="aures")["S"]:.4f} ' + f'(SQAT {p.S_aures:.4f})', end='') + print(f'\n worst |dS| = {worst:.1e} acum') + self.assertLess(worst, 1e-9) + 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)