Big differences in time weighting compared to python-acoustics

[ToGoOrNotToGo]

The results of the time weighting of this library differ completely from those of python-acoustics:
https://github.com/python-acoustics/python-acoustics/blob/master/acoustics/standards/iec_61672_1_2013.py

Is there a reason for this?

[jmrplens]

Hi there, thanks for reporting this.

I have conducted a comprehensive analysis comparing PyOctaveBand with both python-acoustics and the international standard IEC 61672-1:2013 (Electroacoustics - Sound level meters).

Conclusion: PyOctaveBand is implementing the standard correctly. The discrepancies you observed are due to fundamental differences in implementation philosophy and input definitions between the two libraries.

1. Continuous vs. Block Processing: PyOctaveBand provides a continuous time-weighted signal (standard compliant), whereas python-acoustics uses block-based integration (resetting every $\tau$ seconds).
2. Input Units: PyOctaveBand expects raw sound pressure ($Pa$) and handles the squaring internally. python-acoustics expects the "energetic quantity" ($Pa^2$) as input.
3. Stability: PyOctaveBand uses a stable exponential averaging filter. The reference library uses a theoretically unstable filter design that relies on block resetting to function.

---

1. Verification against IEC 61672-1 Standard

The rigorous test for time weighting is the Tone Burst Response (IEC 61672-1, Table 3). I tested both libraries against these targets using a 4 kHz sine burst.

PyOctaveBand Results

Mode (FAST)	Burst Duration	IEC Target (dB)	PyOctaveBand (dB)	Error (dB)	Status
FAST	200 ms	-1.0	-0.98	+0.02	✅ PASS
	50 ms	-4.8	-4.82	-0.02	✅ PASS
	10 ms	-11.1	-11.14	-0.04	✅ PASS
	1 ms	-20.9	-20.99	-0.09	✅ PASS

python-acoustics Results

Note: For these values, I passed the squared signal $x^2$ as required by that library.

Mode (FAST)	Burst Duration	IEC Target (dB)	python-acoustics (dB)	Error (dB)	Status
FAST	200 ms	-1.0	-0.97	+0.03	✅ PASS
	50 ms	-4.8	-3.93	+0.87	⚠️ FAIL
	10 ms	-11.1	-10.90	+0.20	✅ PASS
	1 ms	-20.9	-20.90	+0.00	✅ PASS

PyOctaveBand maintains high precision across all test cases. python-acoustics deviates significantly (>0.8 dB error) for the 50ms burst because its block-based approach ($125ms$ blocks) cannot resolve short transient events accurately depending on their alignment.

---

2. Implementation Differences

A. Continuous vs. Block Output

Standard time weighting is defined as a continuous function of time.

• PyOctaveBand: Returns a continuous envelope.
• python-acoustics: Returns a stepped output (1 value every 125ms).

Step Response Comparison:

B. Theoretical Stability (Detailed)

The standard defines time weighting via the differential equation:

$$ \tau \frac{dy(t)}{dt} + y(t) = x^2(t) $$

This corresponds to a stable Low-Pass Filter with a pole in the Left Half Plane ($s = -1/\tau$).

Interpretation of the Plot:

• PyOctaveBand (Green Dot):

  • Pole location: $s = -8$.
  • Calculation: For FAST weighting, $\tau = 0.125s$. The pole is at $-1/\tau = -1/0.125 = -8$.
  • Physical Meaning: A pole on the negative real axis corresponds to a decaying exponential impulse response ($h(t) \propto e^{-8t}$). This is exactly what "Exponential Time Weighting" means: past events are "forgotten" exponentially. It is an inherently stable system.

• python-acoustics (Red X and O):

  • Pole locations: $s \approx +0.125$ and $s = +1$.
  • Zero location: $s = +1$.
  • Analysis: The zero at $s=1$ theoretically cancels the pole at $s=1$, leaving a dominant pole at $s = +0.125$.
  • Physical Meaning: A pole on the positive real axis (Right Half Plane) corresponds to a growing exponential ($h(t) \propto e^{+0.125t}$). This system is unstable; its output would theoretically grow to infinity over time. The library "fixes" this by forcibly resetting the filter state to zero every 0.125 seconds. This makes it act more like a short-term integrator ($L_{eq, 125ms}$) than a true Time Weighting filter.

---

3. Reproduction Script

You can run the script below to generate these graphs and verify the IEC compliance tables yourself.

Click to expand verification script

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import zpk2tf, bilinear, lfilter
from pyoctaveband import time_weighting

# --- Mimic python-acoustics implementation ---
def acoustics_integrate_mimic(data, sample_frequency, integration_time):
    # This mimics the logic found in the reference library
    b, a = zpk2tf([1.0], [1.0, integration_time], [1.0])
    b, a = bilinear(b, a, fs=sample_frequency)
    n = int(np.floor(integration_time * sample_frequency))
    samples = data.shape[-1]
    # Block processing: Drop tail, reshape, filter, take last value
    data_trimmed = data[..., 0:n * (samples // n)]
    data_reshaped = data_trimmed.reshape(list(data_trimmed.shape[0:-1]) + [-1, n])
    return lfilter(b, a, data_reshaped)[..., n - 1] / integration_time

def check_iec_compliance(fs=48000):
    f0 = 4000
    # IEC 61672-1 Table 3 targets for FAST
    targets = {0.2: -1.0, 0.05: -4.8, 0.01: -11.1, 0.001: -20.9}
    
    print("\n--- IEC 61672-1 Compliance Check (FAST) ---")
    print(f"{ 'Duration':<10} | { 'Target':<8} | { 'PyOctaveBand':<12} | { 'Err(Py)':<8} | { 'Acoustics*':<12} | { 'Err(Ac)':<8}")
    print("-" * 75)
    
    steady_state_val = 0.5
    for dur, target in targets.items():
        total_len = 2.0 
        t = np.linspace(0, total_len, int(fs*total_len), endpoint=False)
        burst_sig = np.zeros_like(t)
        start_idx = int(0.5 * fs)
        end_idx = start_idx + int(dur * fs)
        burst_sig[start_idx:end_idx] = np.sin(2 * np.pi * f0 * t[start_idx:end_idx])
        
        # 1. PyOctaveBand
        y_py = time_weighting(burst_sig, fs, mode="fast")
        val_py = 10 * np.log10(np.max(y_py) / steady_state_val)
        
        # 2. python-acoustics (Squared input)
        y_ac = acoustics_integrate_mimic(burst_sig**2, fs, 0.125)
        val_ac = -99.9
        if len(y_ac) > 0 and np.max(y_ac) > 1e-12:
             val_ac = 10 * np.log10(np.max(y_ac) / steady_state_val)
             
        print(f"{dur*1000:<4.0f} ms    | {target:<8.1f} | {val_py:<12.2f} | {val_py-target:<8.2f} | {val_ac:<12.2f} | {val_ac-target:<8.2f}")
    print("\n*Note: python-acoustics values calculated using squared input signal.")

def plot_comparisons(fs=48000):
    # Step Response
    duration = 0.5
    t = np.linspace(0, duration, int(fs*duration), endpoint=False)
    x = np.sin(2 * np.pi * 1000 * t)
    
    y_py = time_weighting(x, fs, mode="fast")
    y_ac = acoustics_integrate_mimic(x**2, fs, 0.125)
    t_ac = np.arange(len(y_ac)) * 0.125 + 0.125
    
    plt.figure(figsize=(10, 5))
    plt.plot(t, y_py, label='PyOctaveBand (Continuous)')
    plt.plot(t_ac, y_ac, 'o--', label='python-acoustics (Blocks)')
    plt.axhline(0.5, linestyle=':', color='black')
    plt.title("Step Response (Switch On)")
    plt.xlabel("Time [s]")
    plt.ylabel("Mean Square")
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

if __name__ == "__main__":
    check_iec_compliance()
    plot_comparisons()

Conclusions

If your goal is to perform standard-compliant sound level measurements:

• Use PyOctaveBand time_weighting(x) with the raw signal. It provides the correct, stable, and continuous exponential average required by IEC 61672-1.
• If you need to match python-acoustics for legacy reasons, you would need to square your signal first and accept the block-quantized output.

I hope this clarifies the differences!

[mlotinga]

@ToGoOrNotToGo , it's worth noting that python-acoustics is no longer maintained and was forked to an actively maintained repo: https://github.com/Universite-Gustave-Eiffel/acoustic-toolbox — I believe the acoustic-toolbox package should be generally preferred over python-acoustics for this reason.

I think this issue should be raised and linked on that package.

[jmrplens]

Thank you for the information, @mlotinga! I wasn't aware of the acoustic-toolbox fork. I've mostly been focused on bringing my own project back up to speed after a few years of inactivity. I appreciate the heads-up, and whenever I can find some time, I'll definitely look into contributing to that repository.

[ToGoOrNotToGo]

@jmrplens Thank you very much for your excellent answer. I will take a closer look at it in the coming days.

[ToGoOrNotToGo]

The time weighting is an exponential averaging filter, correct? When I compare the results from time_weighting with a direct difference equation implementation, they differ:
y[i]=α⋅x[i]+(1−α)⋅y[i−1]​

Are the results of time_weighting correct when signal.lfilter is used?

[jmrplens]

Hi @ToGoOrNotToGo!

The implementation using signal.lfilter is mathematically identical to the difference equation you proposed. The coefficients used in PyOctaveBand map directly to that recursive formula:

• Difference Equation: $y[i] = \alpha \cdot x^2[i] + (1 - \alpha) \cdot y[i-1]$
• Transfer Function: $H(z) = \frac{\alpha}{1 - (1 - \alpha)z^{-1}}$
• Implementation: b = [α], a = [1, -(1 - α)]

The discrepancies you observed are likely because time_weighting squares the input signal internally (as required by IEC 61672-1). If your manual implementation uses the raw pressure $x$ instead of $x^2$, the results will not match.

Numerical Verification

Running a direct comparison between the lfilter implementation and a manual for loop yields zero numerical difference (see the results below for a 50Hz sine wave at $f_s = 1000$ Hz):

Sample index	PyOctaveBand (lfilter)	Manual Loop (EMA)	Difference
0	0.000000	0.000000	0.0
1	0.000762	0.000762	0.0
2	0.003514	0.003514	0.0
3	0.008708	0.008708	0.0
4	0.015852	0.015852	0.0
...	...	...	...
Max Diff			0.0

Click to expand verification script

import numpy as np
from scipy import signal
from pyoctaveband import time_weighting

def manual_time_weighting(x, fs, mode="fast"):
    # This matches the internal logic of PyOctaveBand
    x_sq = x**2 
    tau = 0.125 if mode == "fast" else 1.0
    alpha = 1 - np.exp(-1 / (fs * tau))
    
    y = np.zeros_like(x_sq)
    curr_y = 0.0
    for i in range(len(x_sq)):
        # EMA Difference Equation: y[i] = alpha * x^2[i] + (1 - alpha) * y[i-1]
        curr_y = alpha * x_sq[i] + (1 - alpha) * curr_y
        y[i] = curr_y
    return y

# Parameters
fs = 1000
t = np.linspace(0, 1, fs)
x = np.sin(2 * np.pi * 50 * t) # 50Hz Sine

# Compare
y_py = time_weighting(x, fs, mode="fast")
y_manual = manual_time_weighting(x, fs, mode="fast")

# Results
diff = np.max(np.abs(y_py - y_manual))
print(f"Mathematical equivalence check: {'PASSED' if diff == 0 else 'FAILED'}")
print(f"Max numerical difference: {diff}")

print("\nFirst 5 samples comparison:")
for i in range(5):
    print(f"Idx {i}: Py={y_py[i]:.6f} | Manual={y_manual[i]:.6f}")

Are you ensuring that the input is squared before applying your manual formula? If you use the raw signal $x$ instead of $x^2$, your results will be in different units (pressure instead of energy), which explains the difference.

[ToGoOrNotToGo]

Thanks for the answer. Now I understand where the difference comes from. An exponential averaging filter is usually initialized with the first value so that it doesn't "swings in" for too long. That's exactly what I had done. In your example:
curr_y = x_sq[0]

I get the same results with curr_y = 0.0. Is it not allowed to initialize the exponential averaging filter when using time-averaging?