diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..b86ef09 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -57,6 +57,11 @@ Modules: * pytta.iso3741: calculations according to the standard; + * pytta.psychoacoustics: + loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel + and Weber 1997), fluctuation strength (Osses et al. 2016) and + tonality (Aures 1985) of a sound pressure signal; + .. toctree:: classes/index utils/index @@ -64,4 +69,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..e7d3c16 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -72,6 +72,11 @@ * pytta.iso3741: calculations according to the standard; + * pytta.psychoacoustics: + loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel + 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. @@ -108,6 +113,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 +136,7 @@ 'generate', 'utils', 'iso3741', + 'psychoacoustics', 'rooms', # Functions diff --git a/pytta/psychoacoustics.py b/pytta/psychoacoustics.py new file mode 100644 index 0000000..e85feb9 --- /dev/null +++ b/pytta/psychoacoustics.py @@ -0,0 +1,1493 @@ +""" +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) + >>> pytta.psychoacoustics.fluctuation_strength_osses2016(signal) + >>> pytta.psychoacoustics.tonality_aures1985(signal) + >>> pytta.psychoacoustics.loudness_level(N) + >>> pytta.psychoacoustics.critical_band(f) + >>> pytta.psychoacoustics.loudness_match(reference, signal) + >>> pytta.psychoacoustics.compare_loudness(signal_a, signal_b) + +* 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). + +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)). + +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 +---------- + +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. + +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. + +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 + +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, +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', + 'fluctuation_strength_osses2016', 'tonality_aures1985', + 'loudness_level', 'loudness_from_level', 'critical_band', + 'loudness_match', 'compare_loudness'] + +_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 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} + + +class _Samples: + # The two attributes of a SignalObj that _mono reads. + def __init__(self, x, fs): + self.timeSignal, self.samplingRate = x, fs + + +def loudness_match(reference, signal, field='free', time_varying=False, + tolerance=0.01): + """ + Gain in dB that gives signal the loudness of reference, as + ita_loudness_adjustment. + + The loudness grows with the gain, so the gain is found by widening a + bracket in steps of 10 dB, as ita_loudness_adjustment does, and solving + on it to 0.001 dB. The loudness is N of loudness_zwicker, or its N5 with + time_varying=True. + + Parameters + ---------- + reference, signal : SignalObj + One channel each, sound pressure in Pa. + field : str, optional + 'free' (default) or 'diffuse'. + time_varying : bool, optional + Compare N5 of the time-varying loudness. Default False. + tolerance : float, optional + Largest relative loudness difference accepted. Default 0.01. + + Returns + ------- + float + Gain in dB to apply to signal. + """ + from scipy.optimize import brentq + method = 'time_varying' if time_varying else 'stationary' + key = 'N5' if time_varying else 'N' + + def loudness(x, fs): + return loudness_zwicker(_Samples(x, fs), field, method)[key] + + target = loudness(*_mono(reference)) + x, fs = _mono(signal) + if fs != _FS: + x, fs = _resample(x, _FS, fs), _FS + if target <= 0 or not np.any(x): + raise ValueError('both signals need a loudness above zero.') + + def excess(g): + return np.log(max(loudness(x * 10 ** (g / 20), fs), 1e-12) / target) + + lo = hi = 0. + f_lo = f_hi = excess(0.) + step = 10. if f_lo < 0 else -10. + while (f_lo < 0) == (f_hi < 0): + if abs(hi) > 200: + raise ValueError('no gain within 200 dB matches the loudness.') + lo, f_lo = hi, f_hi + hi += step + f_hi = excess(hi) + g = brentq(excess, min(lo, hi), max(lo, hi), xtol=1e-3) + if abs(np.expm1(excess(g))) > tolerance: + raise ValueError('the loudness could not be matched within the tolerance.') + return float(g) + + +def compare_loudness(signal_a, signal_b, field='free'): + """ + Ratio of the time-varying specific loudness of two signals, folded to + 1 or more, as ita_loudness_compare_rir. + + Each value is N'_b / N'_a or its inverse, whichever is larger, per 2 ms + frame and 0.1 Bark; it is 1 where both are zero and inf where only one + is, as in ita_loudness_compare_rir. + + Returns + ------- + dict + 'ratio' (frames, 240), 'time' (frames,) and 'bark' (240,). + """ + a = loudness_zwicker(signal_a, field, 'time_varying') + b = loudness_zwicker(signal_b, field, 'time_varying') + if a['specific'].shape != b['specific'].shape: + raise ValueError('both signals need the same duration.') + with np.errstate(divide='ignore', invalid='ignore'): + r = np.abs(b['specific'] / a['specific']) + r = np.where(r < 1, 1 / r, r) + r[np.isnan(r)] = 1 + return {'ratio': r, 'time': a['time'], 'bark': a['bark']} + + +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 + + +_WEIGHTINGS = ('din', 'aures', 'bismarck') + + +def _sharpness(specific, weighting='din'): + """ + 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 + 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, + 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. 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): + ----------------------------- + + * 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; + + * weighting ('din'), (str): + 'din' (Formula (1)), 'aures' (B.1) or 'bismarck' (B.2). + + 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. + + """ + 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'], weighting) + 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 _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 + 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) + 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,)), + (_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] + 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): + """ + 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):]))} + + +# --------------------------------------------------------------------------- +# 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) +# --------------------------------------------------------------------------- + +_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/fluctuation_sqat.mat b/tests/golden/fluctuation_sqat.mat new file mode 100644 index 0000000..964c278 Binary files /dev/null and b/tests/golden/fluctuation_sqat.mat differ diff --git a/tests/golden/loudness_sqat.mat b/tests/golden/loudness_sqat.mat new file mode 100644 index 0000000..3fc7e65 Binary files /dev/null and b/tests/golden/loudness_sqat.mat differ diff --git a/tests/golden/make_fluctuation_sqat.m b/tests/golden/make_fluctuation_sqat.m new file mode 100644 index 0000000..75bffbb --- /dev/null +++ b/tests/golden/make_fluctuation_sqat.m @@ -0,0 +1,94 @@ +function make_fluctuation_sqat(sqatRoot) +% MAKE_FLUCTUATION_SQAT Golden data for tests/test_psychoacoustics.py +% (fluctuation strength). +% +% make_fluctuation_sqat(sqatRoot) runs FluctuationStrength_Osses2016 from +% SQAT (https://github.com/ggrecow/SQAT) on short synthetic signals and +% writes fluctuation_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 paper. 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: +% 'am': A*(1 + m*sin(2*pi*fm*t)).*sin(2*pi*fc*t), +% A = sqrt(2)*2e-5*10^(L/20)/sqrt(1 + m^2/2) (L: rms level); +% 'fm': sqrt(2)*2e-5*10^(L/20)*sin(2*pi*fc*t - (df/fm)*cos(2*pi*fm*t)) +% (frequency fc + df*sin(2*pi*fm*t)); +% +% Fluctuation strength: SQAT with method 1 (time-varying), time_skip 0, +% struct_opt.a0_type = 'fastl2007' (the a0 curve that pytta uses). +% Stored per case (struct array "cases"): the signal parameters, the sum +% and sum of squares of the single precision input, and the outputs FS(t), +% time, specific FS (47 x frames) and FSmean. Also stored: matlab_version and sqat_commit. +% +% Usage: matlab -singleCompThread -batch "make_fluctuation_sqat('/path/to/SQAT')" + +addpath(sqatRoot); +startup_SQAT([strip(sqatRoot, "right", filesep) filesep]); % expects a trailing separator + +fdefs = { +% name kind fc fm m df L fs dur + 'am1k_4hz' 'am' 1000 4 1 0 70 48000 4.0 + 'am1k_16hz' 'am' 1000 16 1 0 70 48000 4.0 + 'am1k_4hz_441' 'am' 1000 4 1 0 60 44100 3.0 + 'fm1k5_4hz' 'fm' 1500 4 0 700 70 48000 4.0 +}; +opt.a0_type = 'fastl2007'; +cases = struct([]); +for c = 1:size(fdefs, 1) + p = cell2struct(fdefs(c, :), {'name','kind','fc','fm','m','df','L','fs','dur'}, 2); + t = (0:round(p.dur * p.fs)-1)' / p.fs; + if strcmp(p.kind, 'am') + A = sqrt(2)*2e-5*10^(p.L/20) / sqrt(1 + p.m^2/2); + x = A * (1 + p.m*sin(2*pi*p.fm*t)) .* sin(2*pi*p.fc*t); + else + x = sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*p.fc*t - (p.df/p.fm)*cos(2*pi*p.fm*t)); + end + x = double(single(x)); + OUT = FluctuationStrength_Osses2016(x, p.fs, 1, 0, 0, opt); + p.sum_x = sum(x); + p.sum_x2 = sum(x.^2); + p.FS = OUT.InstantaneousFluctuationStrength; + p.time = OUT.time; + p.FSspec = OUT.InstantaneousSpecificFluctuationStrength; + p.FSmean = OUT.FSmean; + cases = il_append(cases, p); + fprintf('%-14s FSmean = %.4f vacil (%d frames)\n', p.name, p.FSmean, numel(p.FS)); +end +il_save('fluctuation_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/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/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/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/roughness_sqat.mat b/tests/golden/roughness_sqat.mat new file mode 100644 index 0000000..9b285f5 Binary files /dev/null and b/tests/golden/roughness_sqat.mat differ diff --git a/tests/golden/sharpness_weights_sqat.mat b/tests/golden/sharpness_weights_sqat.mat new file mode 100644 index 0000000..75d7f33 Binary files /dev/null and b/tests/golden/sharpness_weights_sqat.mat differ diff --git a/tests/golden/tonality_sqat.mat b/tests/golden/tonality_sqat.mat new file mode 100644 index 0000000..cbd2974 Binary files /dev/null and b/tests/golden/tonality_sqat.mat differ diff --git a/tests/test_psychoacoustics.py b/tests/test_psychoacoustics.py new file mode 100644 index 0000000..4ee6f60 --- /dev/null +++ b/tests/test_psychoacoustics.py @@ -0,0 +1,565 @@ +""" +Tests for pytta.psychoacoustics. + +Each test pins one rule: a value that defines a unit, a reference value of a +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, +golden/make_roughness_sqat.m, golden/make_fluctuation_sqat.m, +golden/make_tonality_sqat.m and golden/make_sharpness_weights_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 +with the outputs. + +SignalObj inspects the caller's frame (classes/_instanceinfo.py) and fails +when the caller's first local variable is an array, so the objects are built +inline in the tests, or in helpers whose first argument is a scalar. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +import pytta +from pytta import psychoacoustics as pa + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden') +P0 = 2e-5 + + +def golden_cases(name): + d = sio.loadmat(os.path.join(GOLDEN, name), squeeze_me=True, + struct_as_record=False) + return list(np.atleast_1d(d['cases'])) + + +def loudness_signal(p): + """Rebuild the signal of one case of make_loudness_sqat.m.""" + n = int(round(p.dur * p.fs)) + t = np.arange(n) / p.fs + if p.kind in ('tone', 'burst'): + x = np.sqrt(2) * 2e-5 * 10 ** (p.L / 20) * np.sin(2 * np.pi * p.f * t) + if p.kind == 'burst': + x[~((t >= 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 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 + 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 multisine(f1, f2, L, fs=48000, dur=1.0): + """Band noise of make_sharpness_weights_sqat.m: a multisine every 1 Hz + from f1 to f2 with phases 2 pi mod(g k^2, 1), at the rms level L.""" + t = np.arange(int(round(dur * fs))) / fs + k = np.arange(1, int(f2 - f1) + 2) + phases = 2 * np.pi * np.mod(0.6180339887 * k ** 2, 1) + x = np.zeros(t.size) + for kk, ph in zip(k, phases): + x += np.sin(2 * np.pi * (f1 + kk - 1) * t + ph) + return (x * (P0 * 10 ** (L / 20) / np.sqrt(np.mean(x ** 2)))).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 TestScales(unittest.TestCase): + + def test_loudness_level_follows_iso_532_1(self): + # ISO 532-1, Formulae (1) to (3): 1 sone is 40 phon, each doubling + # of N adds 10 phon, and below 1 sone LN = 40 (N + 0.0005)^0.35. + N = np.array([0.01, 0.25, 0.5, 1, 2, 4, 64]) + LN = pa.loudness_level(N) + expected = np.where(N >= 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 TestLoudnessComparison(unittest.TestCase): + + def test_match_returns_the_gain_between_scaled_copies(self): + # A copy 10 dB down needs +10 dB; a different sound gets a gain that + # makes its loudness that of the reference, within the tolerance. + fs = 48000 + t = np.arange(fs) / fs + tone = np.sqrt(2) * P0 * 10 ** (60 / 20) * np.sin(2 * np.pi * 1000 * t) + ref = pytta.SignalObj(tone, 'time', fs) + g = pa.loudness_match(ref, pytta.SignalObj(tone * 10 ** (-10 / 20), 'time', fs)) + noise = multisine(3700, 4400, 50) + g_n = pa.loudness_match(ref, pytta.SignalObj(noise, 'time', fs)) + N_ref = pa.loudness_zwicker(ref)['N'] + N_n = pa.loudness_zwicker(pytta.SignalObj(noise * 10 ** (g_n / 20), 'time', fs))['N'] + print(f'\n scaled copy: {g:.4f} dB; band noise: {g_n:.3f} dB, ' + f'N {N_n:.4f} vs {N_ref:.4f} sone') + self.assertAlmostEqual(g, 10, delta=0.01) + self.assertLess(abs(N_n / N_ref - 1), 0.01) + + def test_compare_is_one_for_identical_sounds_and_symmetric(self): + # The ratio of the specific loudness patterns, folded to >= 1 (1 + # where both are zero), as ita_loudness_compare_rir. + a = multisine(920, 1080, 60, dur=0.5) + b = multisine(920, 1080, 66, dur=0.5) + A, B = (pytta.SignalObj(x, 'time', 48000) for x in (a, b)) + same = pa.compare_loudness(A, A) + ab, ba = pa.compare_loudness(A, B), pa.compare_loudness(B, A) + print(f"\n identical: max |r - 1| {np.abs(same['ratio'] - 1).max():.1e}; " + f"6 dB apart: max r {ab['ratio'].max():.2f}") + self.assertTrue(np.all(same['ratio'] == 1)) + np.testing.assert_allclose(ab['ratio'], ba['ratio']) + self.assertTrue(np.all(ab['ratio'] >= 1)) + self.assertEqual(ab['ratio'].shape, (ab['time'].size, 240)) + + +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_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) + 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 + + +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): + # 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()