[FEATURE] add the sound power level by the direct and comparison methods to pytta.iso3741

pytta/iso3741.py

@@ -5,14 +5,17 @@
 This module does some calculations compliant to ISO 3471 in order to obtain
 the acoustic power of a source.
 
-For now only the calculation of the mean time-average sound pressure level
-among the microphone positions is implemented. In the future we hope the
+Implemented are the mean time-average sound pressure level among the
+microphone positions and the sound power level by the direct and by the
+comparison methods of ISO 3741:1999, 8.4. In the future we hope the
 community keep developing this module in order to provide all the
 calculation needed for measurements acordding to this standard.
 
 Available functions:
 
     >>> pytta.iso3741.Lp_ST()
+    >>> pytta.iso3741.sound_power()
+    >>> pytta.iso3741.sound_power_comparison()
 
 For further information, check the function specific documentation.
 
@@ -140,3 +143,87 @@ def Lp_ST(sigObjList, nthOct, minFreq, maxFreq, IRManualCut=None):
     Lp_ST.error = Sm
     Lp_ST.errorLabel = "Standard deviation"
     return Lp_ST
+
+
+def sound_power(Lp, T, V, S, f, temperature, pressure):
+    """
+    Sound power level by the direct method of ISO 3741:1999, 8.4.1,
+    Equations (11) and (12), normalized to rho c = 400 N s/m3.
+
+        LW = Lp + 10 lg(A/A0) + 4.34 A/S + 10 lg(1 + S c / (8 V f))
+             - 25 lg(427/400 sqrt(273/(273 + theta)) B/B0) - 6 dB
+
+    with A = 55.26 V / (c T), c = 20.05 sqrt(273 + theta) m/s, A0 = 1 m2 and
+    B0 = 1.013e5 Pa. The term 4.34 A/S accounts for air absorption in the
+    room and 10 lg(1 + S c / (8 V f)) is the Waterhouse correction for the
+    energy near the room boundaries.
+
+    Parameters (default), (type):
+    -----------------------------
+
+        * Lp (), (array-like):
+            Average sound pressure level in the room by band, in dB, after
+            the background correction (8.3), e.g. the data of Lp_ST;
+
+        * T (), (array-like):
+            Reverberation time by band with the source installed, in s;
+
+        * V (), (float):
+            Room volume, in m3;
+
+        * S (), (float):
+            Total surface area of the room, in m2;
+
+        * f (), (array-like):
+            Midband frequencies of the bands, in Hz;
+
+        * temperature (), (float):
+            Air temperature in the room, in degC;
+
+        * pressure (), (float):
+            Atmospheric pressure, in Pa.
+
+    Returns:
+    ---------
+
+        numpy.ndarray with LW by band, in dB.
+
+    """
+    c = 20.05 * np.sqrt(273 + temperature)
+    A = 55.26 * V / (c * np.asarray(T, dtype=float))
+    impedance = 427 / 400 * np.sqrt(273 / (273 + temperature)) \
+        * pressure / 1.013e5
+    return np.asarray(Lp, dtype=float) + 10 * np.log10(A) + 4.34 * A / S \
+        + 10 * np.log10(1 + S * c / (8 * V * np.asarray(f, dtype=float))) \
+        - 25 * np.log10(impedance) - 6
+
+
+def sound_power_comparison(Lp, Lp_ref, LW_ref):
+    """
+    Sound power level by the comparison method of ISO 3741:1999, 8.4.2,
+    Equation (13): LW = LWr + (Lp - Lpr).
+
+    Parameters (default), (type):
+    -----------------------------
+
+        * Lp (), (array-like):
+            Average sound pressure level in the room by band with the source
+            under test, in dB (8.3);
+
+        * Lp_ref (), (array-like):
+            Average sound pressure level in the room by band with the
+            reference sound source, in dB (8.4.2.2);
+
+        * LW_ref (), (array-like):
+            Sound power level of the reference sound source by band, in dB,
+            as calibrated for rho c = 400 N s/m3 and not corrected for the
+            conditions of the test.
+
+    Returns:
+    ---------
+
+        numpy.ndarray with LW by band, in dB.
+
+    """
+    return np.asarray(LW_ref, dtype=float) + np.asarray(Lp, dtype=float) \
+        - np.asarray(Lp_ref, dtype=float)

tests/test_iso3741.py

@@ -0,0 +1,59 @@
+"""
+Tests for the sound power level of pytta.iso3741, ISO 3741:1999, 8.4. Each
+test pins one rule.
+"""
+import os
+import unittest
+
+import numpy as np
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import pytta
+
+# A 200 m3 room with 210 m2 of surface, 8 s at 100 Hz and 5 s at 1 000 Hz,
+# at 20 degC and 101 325 Pa.
+ROOM = dict(T=[8.0, 5.0], V=200.0, S=210.0, f=[100.0, 1000.0],
+            temperature=20.0, pressure=101325.0)
+
+
+class TestISO3741SoundPower(unittest.TestCase):
+
+    def test_direct_method_by_hand(self):
+        # Equations (11) and (12), term by term. c = 20.05 sqrt(293)
+        # = 343.20 m/s, A = 55.26 V / (c T) = 4.0253 and 6.4405 m2.
+        # 100 Hz:   10 lg A = 6.0480, 4.34 A/S = 0.0832,
+        #           10 lg(1 + S c / (8 V f)) = 1.6150
+        # 1000 Hz:  8.0892, 0.1331, 0.1914
+        # both:     -25 lg(427/400 sqrt(273/293) 101325/101300) = -0.3281,
+        #           and -6 dB.
+        expected = np.array([1.4182, 2.0856])
+        LW = pytta.iso3741.sound_power([70.0, 70.0], **ROOM)
+        print(f'\nLW - Lp {LW - 70} dB, difference {LW - 70 - expected} dB')
+        np.testing.assert_allclose(LW - 70, expected, atol=1e-4)
+
+    def test_pressure_enters_as_25_lg_rho_c(self):
+        # B appears only in the impedance term, so halving it raises LW by
+        # 25 lg 2 = 7.526 dB in every band.
+        LW = pytta.iso3741.sound_power([70.0, 70.0], **ROOM)
+        low = pytta.iso3741.sound_power([70.0, 70.0],
+                                        **dict(ROOM, pressure=101325.0 / 2))
+        print(f'\nLW(B/2) - LW(B) {low - LW} dB, difference '
+              f'{low - LW - 25 * np.log10(2)} dB')
+        np.testing.assert_allclose(low - LW, 25 * np.log10(2), atol=1e-12)
+
+    def test_comparison_method_agrees_with_the_direct_one(self):
+        # Equation (13). A reference source rated in the same room by the
+        # direct method gives back the direct result of the source under test,
+        # both being normalized to rho c = 400 N s/m3.
+        Lp, Lpr = np.array([72.3, 65.1]), np.array([80.0, 78.4])
+        LWr = pytta.iso3741.sound_power(Lpr, **ROOM)
+        LW = pytta.iso3741.sound_power_comparison(Lp, Lpr, LWr)
+        direct = pytta.iso3741.sound_power(Lp, **ROOM)
+        print(f'\ncomparison - direct {LW - direct} dB')
+        np.testing.assert_allclose(LW, LWr + Lp - Lpr, atol=1e-12)
+        np.testing.assert_allclose(LW, direct, atol=1e-12)
+
+
+if __name__ == '__main__':
+    unittest.main()