Skip to content

Instantly share code, notes, and snippets.

Last active December 13, 2022 04:50
What would you like to do?
A-weighting audio files in Python
#!/usr/bin/env python
# -*- coding: utf-8 -*-
Translated from a MATLAB script (which also includes C-weighting, octave
and one-third-octave digital filters).
Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
Last modification: Aug. 20, 1997, 10:00am.
BSD license
Translated from adsgn.m to Python 2009-07-14
from numpy import pi, polymul
from scipy.signal import bilinear
def A_weighting(fs):
"""Design of an A-weighting filter.
b, a = A_weighting(fs) designs a digital A-weighting filter for
sampling frequency `fs`. Usage: y = scipy.signal.lfilter(b, a, x).
Warning: `fs` should normally be higher than 20 kHz. For example,
fs = 48000 yields a class 1-compliant filter.
[1] IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996.
# Definition of analog A-weighting filter according to IEC/CD 1672.
f1 = 20.598997
f2 = 107.65265
f3 = 737.86223
f4 = 12194.217
A1000 = 1.9997
NUMs = [(2*pi * f4)**2 * (10**(A1000/20)), 0, 0, 0, 0]
DENs = polymul([1, 4*pi * f4, (2*pi * f4)**2],
[1, 4*pi * f1, (2*pi * f1)**2])
DENs = polymul(polymul(DENs, [1, 2*pi * f3]),
[1, 2*pi * f2])
# Use the bilinear transformation to get the digital filter.
# (Octave, MATLAB, and PyLab disagree about Fs vs 1/Fs)
return bilinear(NUMs, DENs, fs)
import sys
from scipy.signal import lfilter
import numpy as np
from A_weighting import A_weighting
import soundfile as sf
except ImportError:
from scikits.audiolab import wavread
# ( can't read 24-bit WAV)
def rms_flat(a): # from matplotlib.mlab
Return the root mean square of all the elements of *a*, flattened out.
return np.sqrt(np.mean(np.absolute(a)**2))
if len(sys.argv) == 2:
filename = sys.argv[1]
filename = 'noise.wav'
x, fs =
except NameError:
x, fs, bits = wavread(filename)
print('Original: {:+.2f} dB'.format(20*np.log10(rms_flat(x))))
b, a = A_weighting(fs)
y = lfilter(b, a, x)
print('A-weighted: {:+.2f} dB'.format(20*np.log10(rms_flat(y))))
Copy link

ThomQu commented Oct 11, 2021


Personnaly, I am trying to build-up a library that wraps up a few key functions for quick audio files analyses. I've picked stuff from librosa, scipy, etc and got your code from some colleague's scripts. That's the reason why I went for that first version.

Thks for the info about the MZTi method, I'll look for a little help around to check what I can do with it.


Copy link


Ok, so maybe I should just make an installable weighting filter library that does nothing else, and then you could use it in yours, etc? I tried combining it with other things like frequency estimation and THD calculation, but maybe those are too unrelated to combine.

Copy link

endolith commented Oct 15, 2021

I was asked on Stack Exchange for my MZTi code, even if crude, so here it is.

I don't think it's correct at very low sample rates, because of aliasing, but should be better at the typical 44.1 or 48 kHz. It's definitely wrong below fs = 2.5 kHz. In between, I think the MZTi zeroes are forcing it into compliance despite being poorly-transformed.

I don't think it makes much practical difference in dB levels vs the BLT method, by the way.

# Poles and zeros for analog filter from specification
z, p, k = ABC_weighting('A')

# Matched Z Transform:
T = 1/fs
z_MZT = np.exp(z*T)  # 4
p_MZT = np.exp(p*T)  # 6

k_MZT = k * ((*pi*1000 - z)/
    *pi*1000 - p)) *
             ( - exp(p*T - 2*pi*1000*T)))/
     - exp(z*T - 2*pi*1000*T)))

# Tweak gain to actual magnitude response TODO: Why is this necessary?
w_MZT, h_MZT = signal.freqz_zpk(z_MZT, p_MZT, k_MZT, [1000/(fs/2)*pi])
k_MZT /= abs(h_MZT)

w = np.array([pi/3, 2*pi/3])
H = abs(signal.freqs_zpk(z, p, k, w*fs)[1] /
        signal.freqz_zpk(z_MZT, p_MZT, k_MZT, w)[1])  # relative magnitude error

# MZTi method of 0, 1/3, 2/3 Nyquist will not work, since response at 0 is 0.
H = np.concatenate(([1], H))

b1 = 0.5 * (H[0] - sqrt(H[0]**2 - 2*H[1]**2 + 2*H[2]**2))
b2 = (3 * (H[0] - b1) -
      sqrt(-3*H[0]**2 + 12*H[1]**2 - 6*H[0]*b1 - 3*b1**2 + 0j)
      ) / 6
b0 = H[0] - b1 - b2
b = np.array([b0, b1, b2])

z_correct = np.roots(b)

# TODO: Is this right?
z_MZTi = np.concatenate((z_MZT, z_correct))
p_MZTi = p_MZT
k_MZTi = k_MZT

# Tweak gain to actual magnitude response again TODO: Why?
w_MZTi, h_MZTi = signal.freqz_zpk(z_MZTi, p_MZTi, k_MZTi, [1000/(fs/2)*pi])
k_MZTi /= abs(h_MZTi)

And the reference is

Copy link

so maybe I should just make an installable weighting filter library that does nothing else

and maybe include bilinear as an option so people can compare with the results of other software?

Copy link

FYI, the MZTi version is accurate down to about 2.5 kHz:

redo 8 kHz
redo 8 kHz closeup

redo 3 kHz
redo 3 kHz closeup

redo 2 kHz
redo 2 kHz closeup

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment