From cbdd2a9e1f05096acce28230b420b239f298f707 Mon Sep 17 00:00:00 2001
From: Sergio Aguirre <djsergioaguirre@gmail.com>
Date: Sat, 31 Oct 2026 22:06:07 +0100
Subject: [PATCH] [FIX] make pytta.iso3741.Lp_ST run and average as ISO 3741
 asks

Lp_ST raised a NameError on the undefined firstChNum and, past it, a
TypeError, since it called fractional_octave_frequencies with minFreq and
maxFreq instead of freqRange. The calibration check now reads the first
channel through channels.mapping.

The average over the microphone positions was the arithmetic mean of the
levels. ISO 3741:1999, Equation (9), asks for the energy mean, while the
standard deviation of Equation (6) is taken around the arithmetic mean,
as before.

Tests check both means with two tones at 80 and 70 dB, and that an
uncalibrated signal is rejected.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 pytta/iso3741.py            | 28 ++++++++--------------
 tests/test_iso3741_lp_st.py | 47 +++++++++++++++++++++++++++++++++++++
 2 files changed, 57 insertions(+), 18 deletions(-)
 create mode 100644 tests/test_iso3741_lp_st.py

diff --git a/pytta/iso3741.py b/pytta/iso3741.py
index ff321e0..ec3bd2a 100644
--- a/pytta/iso3741.py
+++ b/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"
diff --git a/tests/test_iso3741_lp_st.py b/tests/test_iso3741_lp_st.py
new file mode 100644
index 0000000..36f1d5b
--- /dev/null
+++ b/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()