From 2e12da7857b1540d6cf8e38c0447dd6db95241b0 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 22:31:39 +0100 Subject: [PATCH] [FEATURE] add pytta.din45681 with the tonal audibility of DIN 45681 tonal_audibility(signal) finds the tonal components of a noise and their audibility DeltaL = LT - LG - av, DIN 45681:2002: A-weighted spectrum averaged over Hann blocks with lines of 2 Hz, masking level of the critical band by iteration, tone level from the lines within 10 dB of the peak, distinctness by width and flank steepness, and the sum of the tones in one critical band. It is a port of ita_tonality of the ITA-Toolbox, whose BSD licence is reproduced in the module header. Three points differ from the MATLAB code and are documented: the frequency of each FFT line (ita_tonality reads each one 2 Hz high), the absolute scale of the levels (6 dB low there for tones and noise alike, so DeltaL is the same), and the ends of the search for a tone region, which stop at 100 Hz and 8192 Hz. Tests check a tone in white noise against the formula, parity with ita_tonality on golden data from tests/golden/make_din45681_ita.m (DeltaL within 0.005 dB, frequencies 2 Hz below ITA), that noise alone has no tone, and the rejection of invalid input. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/din45681.rst | 7 + docs/pyttadoc/index.rst | 4 + pytta/__init__.py | 5 + pytta/din45681.py | 220 +++++++++++++++++++++++++++++++ tests/golden/din45681_ita.mat | Bin 0 -> 877 bytes tests/golden/make_din45681_ita.m | 31 +++++ tests/test_din45681.py | 110 ++++++++++++++++ 7 files changed, 377 insertions(+) create mode 100644 docs/pyttadoc/din45681.rst create mode 100644 pytta/din45681.py create mode 100644 tests/golden/din45681_ita.mat create mode 100644 tests/golden/make_din45681_ita.m create mode 100644 tests/test_din45681.py diff --git a/docs/pyttadoc/din45681.rst b/docs/pyttadoc/din45681.rst new file mode 100644 index 0000000..02978a2 --- /dev/null +++ b/docs/pyttadoc/din45681.rst @@ -0,0 +1,7 @@ +.. currentmodule:: pytta + +DIN 45681 +========= + +.. automodule:: pytta.din45681 + :members: diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..52d2328 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -57,6 +57,9 @@ Modules: * pytta.iso3741: calculations according to the standard; + * pytta.din45681: + tonal audibility of tones in noise (DIN 45681); + .. toctree:: classes/index utils/index @@ -64,4 +67,5 @@ Modules: functions generate rooms + din45681 properties diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..b77cb10 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -72,6 +72,9 @@ * pytta.iso3741: calculations according to the standard; + * pytta.din45681: + tonal audibility of tones in noise (DIN 45681); + You can find usage examples in the toolbox's examples folder. @@ -108,6 +111,7 @@ from . import utils from . import rooms from . import iso3741 +from . import din45681 from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +134,7 @@ 'generate', 'utils', 'iso3741', + 'din45681', 'rooms', # Functions diff --git a/pytta/din45681.py b/pytta/din45681.py new file mode 100644 index 0000000..48150ff --- /dev/null +++ b/pytta/din45681.py @@ -0,0 +1,220 @@ +""" +PyTTa DIN 45681 tonal audibility: +--------------------------------- + +Detection of tonal components in noise and their audibility, DIN 45681:2002, +as computed by ita_tonality of the ITA-Toolbox. + +Available functions: + + >>> pytta.din45681.tonal_audibility(signal) + +Origin of the code +------------------ + +tonal_audibility is a port to Python of ita_tonality.m of the ITA-Toolbox +(applications/Psychoacoustics), by Daniel Cragg, Institute of Technical +Acoustics, RWTH Aachen University (2010). The method is the same; three +points differ from the MATLAB code and are stated in the docstring: the +frequency of each FFT line (ita_tonality labels line k as (k + 1) df, so +every frequency is 2 Hz high), the absolute scale of the levels (6 dB low +there, for tones and noise alike, so DeltaL does not change), and the ends of +the search for a tone region, which stop at the analysed range. The licence +of the ITA-Toolbox is reproduced below. + +This product includes software developed by the Institute of Technical +Acoustics (RWTH Aachen University). + +References +---------- + +DIN 45681:2002. Akustik - Bestimmung der Tonhaltigkeit von Geraeuschen und +Ermittlung eines Tonzuschlages fuer die Beurteilung von Geraeuschimmissionen. + +Licence of the ported code (ITA-Toolbox, license.txt) +----------------------------------------------------- + + Copyright (c) 2011 + Institute of Technical Acoustics (RWTH Aachen University) + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + 2. 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. + 3. All advertising materials mentioning features or use of this + software must display the following acknowledgment: + + This product includes software developed by the + Institute of Technical Acoustics (RWTH Aachen University). + + 4. Neither the name of the Institute of Technical Acoustics (RWTH + Aachen University), 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 INSTITUTE OF TECHNICAL ACOUSTICS + (RWTH AACHEN UNIVERSITY) ''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 INSTITUTE OF TECHNICAL ACOUSTICS (RWTH AACHEN + UNIVERSITY) 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. +""" + +import numpy as np + +from pytta.classes import SignalObj + +_P0 = 2e-5 +_DF = 2.0 # line spacing, Hz +_FMIN, _FMAX = 100., 8192. +_HANN = 10 * np.log10(1.5) # equivalent noise bandwidth of the Hann window + + +def _a_weighting(f): + # A weighting as written in ita_tonality (12200 Hz pole, 165.5 dB). + return 165.5 + 20 * np.log10(f**4 / (f**2 + 20.6**2) / (f**2 + 12200**2) + / np.sqrt(f**2 + 107.7**2) + / np.sqrt(f**2 + 737.9**2)) + + +def tonal_audibility(signal): + """ + Tonal components of a noise and their audibility, DIN 45681:2002. + + The A-weighted power spectrum is averaged over Hann windowed blocks of + fs / 2 samples (lines of 2 Hz), with a hop of 10000 samples at 44.1 kHz + (the same fraction of the block at other rates). For each line from + 100 Hz to 8192 Hz the masking noise level LS of its critical band, + dfc = 25 + 75 (1 + 1.4 (f/1 kHz)^2)^0.69 Hz, is the energy mean of the + lines not more than 6 dB above LS, iterated until it settles. A + contiguous region of lines 6 dB or more above LS is a tone candidate at + its highest line; its level LT adds the neighbour lines within 10 dB of + the peak. With LG = LS + 10 lg(dfc / df) and the audibility threshold + av = -2 - lg(1 + (f/502)^2.5), the audibility is + + DeltaL = LT - LG - av. + + A candidate is no tone, and gets DeltaL = 0, when it is wider than + 18 (1 + 0.001 f) Hz or when one of its flanks is less steep than 24. + Tones in one critical band are summed into DeltaL_band. + + Differences from ita_tonality: frequencies are those of the FFT lines + (ita_tonality reads each one 2 Hz high); levels are absolute (6 dB low + there for tones and noise, which leaves DeltaL as it is); the search for + a region stops at 100 Hz and 8192 Hz. The tone adjustment K_T of the + standard is not computed. + + Parameters + ---------- + signal : SignalObj + One channel of sound pressure in Pa, at least one block (0.5 s). + + Returns + ------- + dict + 'deltaL': the largest of DeltaL and DeltaL_band (dB), 0 without tones; + 'frequency': the frequency of that tone (Hz), None without tones; + 'tones': one dict per tone with 'frequency', 'bandwidth' (dfc), + 'LS', 'LT', 'LG', 'av', 'deltaL' and 'deltaL_band' (dB, Hz). + + >>> out = pytta.din45681.tonal_audibility(recording) + """ + if not isinstance(signal, SignalObj): + raise TypeError('signal must be a SignalObj.') + if signal.numChannels != 1: + raise ValueError('signal must have one channel.') + fs = signal.samplingRate + N = int(round(fs / _DF)) + df = fs / N + hop = int(round(N * 10000 / 22050)) + x = np.asarray(signal.timeSignal[:, 0], dtype=float) + starts = np.arange(0, x.size - N + 1, hop) + if starts.size == 0: + raise ValueError('signal is shorter than one block of 0.5 s.') + w = np.hanning(N + 2)[1:-1] # MATLAB hanning(N) + blocks = np.stack([x[s:s + N] * w for s in starts]) + X = np.abs(np.fft.rfft(blocks, axis=1))[:, :N // 2] + power = np.mean(2 * X**2 / np.sum(w)**2, axis=0) / _P0**2 + freq = np.arange(N // 2) * df + with np.errstate(divide='ignore'): + level = 10 * np.log10(power) + level[1:] += _a_weighting(freq[1:]) + level[0] = -np.inf + + first, last = int(round(_FMIN / df)), int(round(_FMAX / df)) + idx = np.arange(first, last + 1) + dfc = 25 + 75 * (1 + 1.4 * (freq / 1000)**2)**0.69 + lo = np.clip(np.searchsorted(freq, freq - dfc / 2), 1, last) + hi = np.clip(np.searchsorted(freq, freq + dfc / 2) - 1, None, last) + + LS = np.full(freq.size, np.inf) + for i in idx: # masking noise level + band = level[lo[i]:hi[i] + 1] + ls, previous, lines = 1000., 0., 1000 + while lines > 9 and previous != ls: + keep = band <= ls + 6 + previous, lines = ls, int(np.sum(keep)) + if lines: + ls = 10 * np.log10(np.mean(10**(band[keep] / 10))) - _HANN + LS[i] = previous if lines < 10 else ls + + above = level >= LS + 6 + deltaL, LT, LG, av = {}, {}, {}, {} + i = first + while i <= last: # tone candidates + if not above[i]: + i += 1 + continue + end = i + while end + 1 <= last and above[end + 1]: + end += 1 + k = i + int(np.argmax(level[i:end + 1])) + ref = LS[k] + 6 + tone, nlines = 10**(level[k] / 10), 1 + j = k - 1 + while level[j] >= ref and level[j] >= level[k] - 10: + tone, nlines, j = tone + 10**(level[j] / 10), nlines + 1, j - 1 + lower = ((level[k] - level[j]) * freq[k] / np.sqrt(2) / (freq[k] - freq[j]) + if j < k - 1 else 100) + j = k + 1 + while level[j] >= ref and level[j] >= level[k] - 10: + tone, nlines, j = tone + 10**(level[j] / 10), nlines + 1, j + 1 + upper = ((level[k] - level[j]) * freq[k] / np.sqrt(2) / (freq[j] - freq[k]) + if j > k + 1 else 100) + LT[k] = 10 * np.log10(tone) - (_HANN if nlines > 1 else 0) + LG[k] = LS[k] + 10 * np.log10(dfc[k] / df) + av[k] = -2 - np.log10(1 + (freq[k] / 502)**2.5) + deltaL[k] = LT[k] - LG[k] - av[k] + if (nlines * df > 18 * (1 + 0.001 * freq[k]) or lower < 24 + or upper < 24): + deltaL[k] = 0. + i = end + 1 + + tones = [] + for k in sorted(k for k in deltaL if deltaL[k] > 0): + others = [m for m in deltaL if m != k and deltaL[m] > 0 + and lo[k] <= m <= hi[k]] + band = (10 * np.log10(sum(10**(LT[m] / 10) for m in [k] + others)) + - LG[k] - av[k]) if others else deltaL[k] + tones.append(dict(frequency=freq[k], bandwidth=dfc[k], LS=LS[k], + LT=LT[k], LG=LG[k], av=av[k], deltaL=deltaL[k], + deltaL_band=band)) + if not tones: + return dict(deltaL=0, frequency=None, tones=[]) + best = max(tones, key=lambda t: max(t['deltaL'], t['deltaL_band'])) + return dict(deltaL=max(best['deltaL'], best['deltaL_band']), + frequency=best['frequency'], tones=tones) diff --git a/tests/golden/din45681_ita.mat b/tests/golden/din45681_ita.mat new file mode 100644 index 0000000000000000000000000000000000000000..c9f05890cd129c77c7ab03c906adbc608698aba4 GIT binary patch literal 877 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X01nwjV*I2WZRmZYXABu-G3Rk|Lc)!MM-Loe zTVmNTqmz*-VquM&LR#Y`g^f>+Ji3x(FuS#pTT71hNM+~Fl>&2*u3V`p5EfULR+v>2 zm}lthStl1E&g|@$6`931BQjd+F-2wFC8izLr z8Vc08#mYZ9W(o*2eBu?|WV=el=(EUyzQztog_(!jj>@>H6wheTnQ>E0<;CNcuJG

d%4Pf$C**Gy!`bVX;?6Wf#D*A;!WkH5cL>EGO)%$EOZS-<8pGuFoN ztoX*}aFd7eFWyHqyS*vgV6f=fp!_mIb{05Pkfj;-{=F&416F|8}&ui?98<qFl_oQEAe5udCt?>hN7p;)7lz6^R0KBEM70j`mWR> zeTU`##`}#if4ASS?XEh%#4D`o_EF!&D(%R33u5cl(vS8>v0pp-)q2hOS5|jSxT211 zxN!F2uVssWx=uXDXScg`Zocm$dzio2f&NlE1@hOQNuO$13v&6qiXEa9H4XJ!8~0SC zBrfpZmF(2XE|#!Y>!{NzJ2eFpbM^as*)N$(*7^MOFWYMsyfAc&%GsVpRbQFpY8Cb` zTE3e#z-&bv>t$Bq2Vb*8V|U%EI~Ky;zo}e)``*h%J8I|Ezgoz8i|-7ZQCrztIlj-| z#lOxv^QquN^o@n*o<4l=;Z|f?bfER3f3YjJna%!gd%JP|xz>Nu%ssL>&oquMF_LzE zlPJ(|`Cw62@RQa(4@`7qrn!`*Mb(>?GB6m{|Ih@bUUeXbq+SDtBr}C)4VM%Sa6L6X YWGL($Wys9TJ)tb}@Bsz}c~2%20GzsU^8f$< literal 0 HcmV?d00001 diff --git a/tests/golden/make_din45681_ita.m b/tests/golden/make_din45681_ita.m new file mode 100644 index 0000000..2a18c4a --- /dev/null +++ b/tests/golden/make_din45681_ita.m @@ -0,0 +1,31 @@ +% Golden data for tests/test_din45681.py: ita_tonality (ITA-Toolbox, +% applications/Psychoacoustics) on tones in white noise at 44.1 kHz. +% The noise is uniform, from rand with the Mersenne twister, so that +% numpy.random.RandomState(seed).random_sample rebuilds the same samples; +% the sum of each signal is stored to check that. +% Run from tests/golden with the ITA-Toolbox on the path. +fs = 44100; t = (0:3*fs-1)'/fs; p0 = 2e-5; +cases = struct('name', {}, 'seed', {}, 'noise_rms', {}, 'tones', {}, 'levels', {}); +cases(1) = struct('name', '1 kHz in noise', 'seed', 1, 'noise_rms', 0.02, 'tones', 1000, 'levels', 60); +cases(2) = struct('name', '400 Hz in noise', 'seed', 2, 'noise_rms', 0.02, 'tones', 400, 'levels', 55); +cases(3) = struct('name', '3 kHz in noise', 'seed', 3, 'noise_rms', 0.02, 'tones', 3000, 'levels', 50); +cases(4) = struct('name', 'two tones in one band', 'seed', 4, 'noise_rms', 0.02, 'tones', [1000 1050], 'levels', [55 55]); +for k = 1:numel(cases) + c = cases(k); + rng(c.seed, 'twister'); + n = rand(numel(t), 1) - 0.5; + x = c.noise_rms * n / sqrt(1/12); + for m = 1:numel(c.tones) + x = x + sqrt(2) * p0 * 10^(c.levels(m)/20) * sin(2*pi*c.tones(m)*t); + end + x = double(single(x)); + a = itaAudio(x, fs, 'time'); a.channelUnits = {'Pa'}; + [dl, f, info] = ita_tonality(a); + close all + cases(k).sum = sum(x); + cases(k).deltaL = dl; + cases(k).frequency = f; + cases(k).info = info; + fprintf('%s: deltaL %.3f dB at %g Hz, %d tones\n', c.name, dl, f, size(info, 1)); +end +save('din45681_ita.mat', 'cases', 'fs'); diff --git a/tests/test_din45681.py b/tests/test_din45681.py new file mode 100644 index 0000000..c4a921f --- /dev/null +++ b/tests/test_din45681.py @@ -0,0 +1,110 @@ +""" +Tests for pytta.din45681, the tonal audibility of DIN 45681:2002 as +computed by ita_tonality of the ITA-Toolbox. Each test pins one rule. + +The ITA golden data come from golden/make_din45681_ita.m. The noise there +is uniform noise from MATLAB's Mersenne twister, which +numpy.random.RandomState rebuilds; the stored sums check the samples. + +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 after a scalar local. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import pytta +from pytta import din45681 + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden') +P0 = 2e-5 + + +def golden_signal(case, fs): + """Rebuild the signal of one case of make_din45681_ita.m.""" + t = np.arange(3 * fs) / fs + noise = np.random.RandomState(int(case.seed)).random_sample(t.size) - 0.5 + x = case.noise_rms * noise / np.sqrt(1 / 12) + for f, L in zip(np.atleast_1d(case.tones), np.atleast_1d(case.levels)): + x = x + np.sqrt(2) * P0 * 10**(L / 20) * np.sin(2 * np.pi * f * t) + return x.astype(np.float32).astype(float) + + +def a_weighting(f): + return 20 * np.log10(12194**2 * f**4 / ((f**2 + 20.6**2) * np.sqrt( + (f**2 + 107.7**2) * (f**2 + 737.9**2)) * (f**2 + 12194**2))) + 2.0 + + +class TestTonalAudibility(unittest.TestCase): + + def test_tone_in_white_noise_against_the_formula(self): + # DeltaL = LT - LG - av, LG = LS + 10 lg(dfc / df), with LS the + # level of white noise in one 2 Hz line and av = -2 - lg(1 + + # (f/502)^2.5). At 1 kHz the A weighting is 0 dB, so the tone keeps + # its level; the noise lines around it keep theirs within 0.1 dB. + # LS is estimated from about 80 lines of 19 blocks: over 8 noise + # seeds it scatters from -0.22 to +0.34 dB, hence 0.5 dB. + label = 'tone in noise' + fs, f, L, sigma = 48000, 1000.0, 60.0, 0.02 + t = np.arange(4 * fs) / fs + x = (sigma * np.random.RandomState(7).standard_normal(t.size) + + np.sqrt(2) * P0 * 10**(L / 20) * np.sin(2 * np.pi * f * t)) + out = din45681.tonal_audibility(pytta.SignalObj(x, 'time', fs)) + N = fs // 2 # lines of 2 Hz; one line holds 2 sigma^2 / N of the noise + LS = 10 * np.log10(2 * sigma**2 / N / P0**2) + a_weighting(f) + dfc = 25 + 75 * (1 + 1.4 * (f / 1000)**2)**0.69 + expected = L - (LS + 10 * np.log10(dfc / 2)) + 2 + np.log10(1 + (f / 502)**2.5) + print(f'\n DeltaL = {out["deltaL"]:.3f} dB at {out["frequency"]:.0f} Hz, ' + f'expected {expected:.3f} dB, diff {out["deltaL"] - expected:+.3f}') + self.assertAlmostEqual(out['frequency'], f, delta=2) + self.assertAlmostEqual(out['deltaL'], expected, delta=0.5) + + def test_matches_ita_tonality(self): + # Same samples in both. ita_tonality reads every frequency 2 Hz high + # (its frequency vector starts at 2 Hz while the first FFT line is + # 0 Hz), so its frequencies are compared less 2 Hz. The rest of the + # method is the same; the A weighting read 2 Hz off shifts DeltaL by + # far less than 0.05 dB. + d = sio.loadmat(os.path.join(GOLDEN, 'din45681_ita.mat'), + squeeze_me=True, struct_as_record=False) + fs = int(d['fs']) + for case in np.atleast_1d(d['cases']): + label = case.name + x = golden_signal(case, fs) + self.assertAlmostEqual(np.sum(x), case.sum, delta=1e-9) + out = din45681.tonal_audibility(pytta.SignalObj(x, 'time', fs)) + info = np.atleast_2d(case.info) + print(f'\n {case.name}: DeltaL {out["deltaL"]:.3f} vs {case.deltaL:.3f} dB, ' + f'{out["frequency"]:.0f} vs {case.frequency - 2:.0f} Hz, ' + f'{len(out["tones"])} vs {info.shape[0]} tones', end='') + self.assertAlmostEqual(out['deltaL'], case.deltaL, delta=0.05) + self.assertAlmostEqual(out['frequency'], case.frequency - 2, delta=1e-9) + self.assertEqual(len(out['tones']), info.shape[0]) + + def test_noise_has_no_tone(self): + label = 'noise' + fs = 48000 + x = 0.02 * np.random.RandomState(3).standard_normal(4 * fs) + out = din45681.tonal_audibility(pytta.SignalObj(x, 'time', fs)) + print(f'\n {len(out["tones"])} tones, DeltaL {out["deltaL"]}') + self.assertEqual(out['tones'], []) + self.assertEqual(out['deltaL'], 0) + + def test_rejects_invalid_input(self): + label = 'invalid' + fs = 48000 + with self.assertRaises(ValueError): + din45681.tonal_audibility(pytta.SignalObj(np.zeros((4 * fs, 2)), 'time', fs)) + with self.assertRaises(ValueError): + din45681.tonal_audibility(pytta.SignalObj(np.zeros(fs // 4), 'time', fs)) + with self.assertRaises(TypeError): + din45681.tonal_audibility(np.zeros(4 * fs)) + + +if __name__ == '__main__': + unittest.main()