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 0000000..c9f0589 Binary files /dev/null and b/tests/golden/din45681_ita.mat differ 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()