[FEATURE] add pytta.iso10534 for impedance tube measurements

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

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))

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()