diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..31433fb 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -108,6 +108,7 @@ from . import utils from . import rooms from . import iso3741 +from . import iso10534 from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +131,7 @@ 'generate', 'utils', 'iso3741', + 'iso10534', 'rooms', # Functions diff --git a/pytta/iso10534.py b/pytta/iso10534.py new file mode 100644 index 0000000..aaa7fbe --- /dev/null +++ b/pytta/iso10534.py @@ -0,0 +1,68 @@ +""" +PyTTa ISO 10534-2 impedance tube: +-------------------------------- +Normal-incidence reflection factor, absorption coefficient and impedance +of a sample in an impedance tube, by the two-microphone transfer function +method of ISO 10534-2:1998. + + >>> f, H12 = pytta.iso10534.transfer_function(mic1, mic2) + >>> r = pytta.iso10534.reflection_factor(H12, f, x1, s, c0) + >>> alpha = pytta.iso10534.absorption(r) + +Microphone 1 is the one further from the sample, at x1 from it; microphone +2 is s closer. The time convention is that of numpy.fft, exp(+j w t), the +same as the standard's. +""" +import numpy as np + + +def speed_of_sound(temperature): + """Speed of sound in m/s at `temperature` in degrees Celsius, + Equation (5): c0 = 343.2 sqrt(T / 293), T in kelvin.""" + return 343.2 * np.sqrt((temperature + 273.15) / 293) + + +def frequency_range(d, s, c0): + """ + Working frequency range (f_low, f_up) in Hz of a circular tube of inner + diameter d with microphones s apart, both in m: f_up from Equations (2) + and (4) (f_up d < 0.58 c0, f_up s < 0.45 c0), and f_low where s is 5 % + of the wavelength, the general guide of 4.2. + """ + return 0.05 * c0 / s, min(0.58 * c0 / d, 0.45 * c0 / s) + + +def transfer_function(mic1, mic2): + """Frequencies and transfer function H12 = P2 / P1 from the recordings + or impulse responses of the two microphone positions (first channel of + each SignalObj).""" + p1 = np.asarray(mic1.timeSignal[:, 0], dtype=float) + p2 = np.asarray(mic2.timeSignal[:, 0], dtype=float) + return np.fft.rfftfreq(p1.size, 1 / mic1.samplingRate)[1:], \ + (np.fft.rfft(p2) / np.fft.rfft(p1))[1:] + + +def reflection_factor(H12, f, x1, s, c0, calibration=None): + """ + Normal-incidence reflection factor, Equation (17): + r = (H12 - HI) / (HR - H12) exp(2 j k x1), with HI = exp(-j k s) and + HR = exp(j k s) (D.5, D.6). `calibration` is the factor Hc of 7.5.2, + which divides H12, if the microphones were not interchanged. + """ + H12 = np.asarray(H12) if calibration is None \ + else np.asarray(H12) / calibration + k = 2 * np.pi * np.asarray(f) / c0 + return (H12 - np.exp(-1j * k * s)) / (np.exp(1j * k * s) - H12) \ + * np.exp(2j * k * x1) + + +def absorption(r): + """Normal-incidence absorption coefficient, Equation (18): + 1 - |r|^2.""" + return 1 - np.abs(r) ** 2 + + +def impedance_ratio(r): + """Specific acoustic impedance ratio Z / (rho c0), Equation (19): + (1 + r) / (1 - r).""" + return (1 + np.asarray(r)) / (1 - np.asarray(r)) diff --git a/tests/test_iso10534.py b/tests/test_iso10534.py new file mode 100644 index 0000000..d4f32ab --- /dev/null +++ b/tests/test_iso10534.py @@ -0,0 +1,67 @@ +""" +Tests for pytta.iso10534, the two-microphone transfer function method of +ISO 10534-2 for impedance tubes. The transfer functions are synthesized +with Equation (D.7) for known reflection factors, so the expected values +are exact. Each test pins one rule. +""" +import os +import unittest + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import numpy as np + +import pytta +from pytta import iso10534 + +F = np.linspace(100, 1800, 200) +C = iso10534.speed_of_sound(20.0) +X1, S = 0.10, 0.05 # m: sample to far microphone, spacing + + +def h12(r): + """Equation (D.7): H12 = p2 / p1 for a reflection factor r at x = 0.""" + k = 2 * np.pi * F / C + x2 = X1 - S + return (np.exp(1j * k * x2) + r * np.exp(-1j * k * x2)) \ + / (np.exp(1j * k * X1) + r * np.exp(-1j * k * X1)) + + +class TestISO10534(unittest.TestCase): + + def test_speed_of_sound(self): + # Equation (5): c0 = 343.2 sqrt(T / 293), T in kelvin. + self.assertAlmostEqual(iso10534.speed_of_sound(19.85), 343.2, 6) + + def test_known_reflection_factors_come_back(self): + # Rigid termination (r = 1, alpha = 0), anechoic (r = 0, alpha = 1, + # Z = rho c) and a complex r with frequency-dependent phase. + cases = (np.ones(F.size), np.zeros(F.size), + 0.5 * np.exp(1j * (np.pi / 4 + F / 1000))) + for r in cases: + got = iso10534.reflection_factor(h12(r), F, X1, S, C) + self.assertTrue(np.allclose(got, r, atol=1e-9)) + alpha = iso10534.absorption(got) + self.assertTrue(np.allclose(alpha, 1 - np.abs(r) ** 2, atol=1e-9)) + z = iso10534.impedance_ratio(np.zeros(F.size)) + self.assertTrue(np.allclose(z, 1)) + + def test_transfer_function_from_two_responses(self): + # p2 is p1 delayed by d samples and halved: H12 = 0.5 exp(-j w d). + fs, d = 48000, 7 + p1 = pytta.SignalObj(np.eye(1, 4096, 100)[0], 'time', fs) + p2 = pytta.SignalObj(0.5 * np.eye(1, 4096, 100 + d)[0], 'time', fs) + f, H = iso10534.transfer_function(p1, p2) + self.assertTrue(np.allclose(H, 0.5 * np.exp(-2j * np.pi * f * d / fs))) + + def test_working_frequency_range(self): + # Equations (2) and (4) and the 5 % rule: a 100 mm circular tube + # with microphones 50 mm apart. + c = 343.2 + low, high = iso10534.frequency_range(0.1, 0.05, c) + self.assertAlmostEqual(high, min(0.58 * c / 0.1, 0.45 * c / 0.05)) + self.assertAlmostEqual(low, 0.05 * c / 0.05) + + +if __name__ == '__main__': + unittest.main()