[FEATURE] add pytta.iso11654 with the weighted sound absorption coefficient

docs/pyttadoc/index.rst

@@ -64,4 +64,5 @@ Modules:
    functions
    generate
    rooms
+   iso11654
    properties

docs/pyttadoc/iso11654.rst

@@ -0,0 +1,7 @@
+.. currentmodule:: pytta
+
+ISO 11654
+=========
+
+.. automodule:: pytta.iso11654
+   :members:

pytta/__init__.py

@@ -72,6 +72,9 @@
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.iso11654:
+        rating of sound absorption (ISO 11654);
+
 
 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 iso11654
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +134,7 @@
            'generate',
            'utils',
            'iso3741',
+           'iso11654',
            'rooms',
 
            # Functions

pytta/iso11654.py

@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+"""
+Rating of sound absorption, ISO 11654:1997: the practical sound absorption
+coefficients of the octave bands from one-third octave values measured by
+ISO 354, the weighted sound absorption coefficient alpha_w and the shape
+indicators L, M and H.
+
+>>> ap = pytta.iso11654.practical_absorption(alpha_third_octaves)
+>>> alpha_w, shape = pytta.iso11654.weighted_absorption(ap)
+
+"""
+
+import numpy as np
+
+__all__ = ['OCTAVES', 'REFERENCE', 'practical_absorption', 'weighted_absorption',
+           'rating']
+
+OCTAVES = np.array([125, 250, 500, 1000, 2000, 4000])
+# Reference curve, 250 Hz to 4 kHz (Figure 1).
+REFERENCE = np.array([0.80, 1.00, 1.00, 1.00, 0.90])
+
+
+def practical_absorption(alpha):
+    """
+    Practical sound absorption coefficients alpha_p of the octave bands, one
+    per three consecutive one-third octave values of alpha (4.1): their
+    arithmetic mean to the second decimal, rounded in steps of 0.05 and
+    limited to 1.00. Give the 18 one-third octave bands from 100 Hz to 5 kHz
+    for the six octaves from 125 Hz to 4 kHz.
+    """
+    a = np.asarray(alpha, dtype=float)
+    if a.size % 3:
+        raise ValueError('Give three one-third octave values per octave.')
+    hundredths = np.round(a.reshape(-1, 3).mean(axis=1) * 100)
+    return np.minimum(np.round(hundredths / 5) * 5 / 100, 1.0)
+
+
+def weighted_absorption(alpha_p):
+    """
+    Weighted sound absorption coefficient alpha_w and the shape indicators
+    (4.2 and 4.3) from the practical coefficients of the octaves 250 Hz to
+    4 kHz (5 values, or 6 starting at 125 Hz). The reference curve moves in
+    steps of 0.05 towards the measured values until the sum of the
+    unfavourable deviations, where the measured value is below the curve,
+    is 0.10 or less; alpha_w is the shifted curve at 500 Hz. A practical
+    coefficient 0.25 or more above the shifted curve adds L (250 Hz), M (500
+    Hz or 1 kHz) or H (2 kHz or 4 kHz). Returns (alpha_w, indicators).
+    """
+    a = np.asarray(alpha_p, dtype=float)
+    if a.size == 6:
+        a = a[1:]
+    if a.size != 5:
+        raise ValueError('Give the octaves 250 Hz to 4 kHz (or 125 Hz to 4 kHz).')
+    # Work in hundredths, so that the steps of 0.05 are exact.
+    ap = np.round(a * 100).astype(int)
+    ref = np.round(REFERENCE * 100).astype(int)
+    shift = -100
+    while np.sum(np.maximum(ref + shift + 5 - ap, 0)) <= 10:
+        shift += 5
+    curve = ref + shift
+    excess = ap - curve >= 25
+    shape = ''.join(letter for letter, bands in (('L', [0]), ('M', [1, 2]), ('H', [3, 4]))
+                    if excess[bands].any())
+    return curve[1] / 100, shape
+
+
+def rating(alpha_w, shape):
+    """The rating as ISO 11654 writes it, e.g. 0.60(M); two decimals."""
+    return f'{alpha_w:.2f}' + (f'({shape})' if shape else '')

tests/test_iso11654.py

@@ -0,0 +1,37 @@
+"""
+Tests for pytta.iso11654 with the rules and the examples of ISO 11654:1997.
+"""
+import unittest
+
+import numpy as np
+
+from pytta import iso11654
+
+
+class TestISO11654(unittest.TestCase):
+
+    def test_rounding_of_the_practical_coefficients(self):
+        # 4.1: x,y2 -> x,y0, x,y3 -> x,y5, x,y7 -> x,y5, x,y8 -> x,y + 0,1,
+        # 0,92 -> 0,90, and at most 1,00.
+        means = [0.42, 0.43, 0.47, 0.48, 0.92, 1.07]
+        third = np.repeat(means, 3)
+        np.testing.assert_allclose(iso11654.practical_absorption(third),
+                                   [0.40, 0.45, 0.45, 0.50, 0.90, 1.00])
+
+    def test_examples_of_annex_a(self):
+        # Figure A.1: alpha_w = 0,60 with an unfavourable deviation at 250
+        # Hz; Figure A.2: the same with 1,00 at 500 Hz gives 0,60(M).
+        a1 = [0.35, 0.70, 0.65, 0.60, 0.55]
+        a2 = [0.35, 1.00, 0.65, 0.60, 0.55]
+        self.assertEqual(iso11654.weighted_absorption(a1), (0.60, ''))
+        self.assertEqual(iso11654.weighted_absorption(a2), (0.60, 'M'))
+        self.assertEqual(iso11654.rating(*iso11654.weighted_absorption(a2)), '0.60(M)')
+
+    def test_full_absorber_and_shape_indicators(self):
+        # All ones rate 1.00; a resonator absorbing at 250 Hz only gets L.
+        self.assertEqual(iso11654.weighted_absorption([1.0] * 6), (1.00, ''))
+        self.assertEqual(iso11654.weighted_absorption([0.9, 0.2, 0.2, 0.2, 0.2]), (0.20, 'L'))
+
+
+if __name__ == '__main__':
+    unittest.main()