diff --git a/pytta/iso3741.py b/pytta/iso3741.py index ff321e0..d81875e 100644 --- a/pytta/iso3741.py +++ b/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) diff --git a/tests/test_iso3741.py b/tests/test_iso3741.py new file mode 100644 index 0000000..2a8aa41 --- /dev/null +++ b/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()