[FIX] make pytta.iso3741.Lp_ST run and average as ISO 3741 asks

pytta/iso3741.py

@@ -94,7 +94,7 @@ def Lp_ST(sigObjList, nthOct, minFreq, maxFreq, IRManualCut=None):
     Leqs = []
 
     for idx, sigObj in enumerate(sigObjList):
-        if not sigObj.channels[firstChNum].calibCheck:
+        if not sigObj.channels[sigObj.channels.mapping[0]].calibCheck:
             raise ValueError("SignalObj {} must be calibrated.".format(idx+1))
         # Cutting the IR
         if IRManualCut is not None:
@@ -106,8 +106,7 @@ def Lp_ST(sigObjList, nthOct, minFreq, maxFreq, IRManualCut=None):
         hSignal = _filter(signal=hSignal, nthOct=nthOct, minFreq=minFreq,
                         maxFreq=maxFreq)
         bands = FOF(nthOct=nthOct,
-                    minFreq=minFreq,
-                    maxFreq=maxFreq)[:,1]
+                    freqRange=[minFreq, maxFreq])[:,1]
         Leq = []
         for chIndex in range(hSignal.numChannels):
             Leq.append(
@@ -119,23 +118,16 @@ def Lp_ST(sigObjList, nthOct, minFreq, maxFreq, IRManualCut=None):
                            comment='Leq')
         Leqs.append(Leq)
 
-    Leq = 0
-    for L in Leqs:
-        Leq =  L + Leq
-    Lp_ST = Leq / len(sigObjList)
-    Lp_ST.anType = 'mixed'
+    # Equation (9): energy mean over the microphone positions; Equation (6):
+    # standard deviation around the arithmetic mean of the levels.
+    data = np.vstack([an.data for an in Leqs])
+    Lp_ST = Analysis(anType='mixed', nthOct=nthOct,
+                     minBand=Leqs[0].minBand, maxBand=Leqs[0].maxBand,
+                     data=10*np.log10(np.mean(10**(data/10), axis=0)),
+                     comment='Lp(ST)')
     Lp_ST.unit = 'dB'
     Lp_ST.creation_name = creation_name
-
-    # Statistics for Lp_ST
-    data = np.vstack([an.data for an in Leqs])
-    Sm = []
-    for bandIdx in range(data.shape[1]):
-        summing = 0
-        for idx in range(data.shape[0]):
-            summing += \
-            (data[idx, bandIdx] - Lp_ST.data[bandIdx])**2 / (data.shape[0] - 1)
-        Sm.append(summing**(1/2))
+    Sm = list(np.std(data, axis=0, ddof=1))
 
     Lp_ST.error = Sm
     Lp_ST.errorLabel = "Standard deviation"

tests/test_iso3741_lp_st.py

@@ -0,0 +1,47 @@
+"""
+Tests for pytta.iso3741.Lp_ST, the average sound pressure level over the
+microphone positions of ISO 3741:1999, 8.3 and 8.1.5. Each test pins one
+rule.
+"""
+import os
+import unittest
+
+import numpy as np
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import pytta
+
+FS = 48000
+
+
+def tone(level, cal=True):
+    """1 kHz tone of `level` dB re 20 uPa, flagged as calibrated."""
+    label = 'tone'
+    t = np.arange(2 * FS) / FS
+    amplitude = np.sqrt(2) * 2e-5 * 10**(level / 20)
+    sig = pytta.SignalObj(amplitude * np.sin(2 * np.pi * 1000 * t), 'time', FS)
+    sig.channels[sig.channels.mapping[0]].calibCheck = cal
+    return sig
+
+
+class TestLpST(unittest.TestCase):
+
+    def test_energy_mean_over_positions_and_arithmetic_deviation(self):
+        # Equation (9): Lp = 10 lg(mean 10^(0.1 Lpi)), so 80 and 70 dB give
+        # 77.40 dB (the arithmetic mean would be 75 dB). Equation (6): sM is
+        # the deviation around the arithmetic mean, 7.07 dB.
+        result = pytta.iso3741.Lp_ST([tone(80), tone(70)], 3, 1000, 1000)
+        expected = 10 * np.log10((1e8 + 1e7) / 2)
+        print(f'\n  Lp = {result.data[0]:.3f} dB (expected {expected:.3f}), '
+              f'sM = {result.error[0]:.3f} dB (expected {np.sqrt(50):.3f})')
+        self.assertAlmostEqual(result.data[0], expected, delta=0.1)
+        self.assertAlmostEqual(result.error[0], np.sqrt(50), delta=0.1)
+
+    def test_rejects_uncalibrated_signal(self):
+        with self.assertRaises(ValueError):
+            pytta.iso3741.Lp_ST([tone(80), tone(70, cal=False)], 3, 1000, 1000)
+
+
+if __name__ == '__main__':
+    unittest.main()