Skip to content

Instantly share code, notes, and snippets.

@endolith
Last active December 13, 2022 04:50
Embed
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)
couvreur@thor.fpms.ac.be
Last modification: Aug. 20, 1997, 10:00am.
BSD license
http://www.mathworks.com/matlabcentral/fileexchange/69
Translated from adsgn.m to Python 2009-07-14 endolith@gmail.com
"""
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.
References:
[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
try:
import soundfile as sf
except ImportError:
from scikits.audiolab import wavread
# (scipy.io.wavfile 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]
else:
filename = 'noise.wav'
try:
x, fs = sf.read(filename)
except NameError:
x, fs, bits = wavread(filename)
print(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))))
@ThomQu
Copy link

ThomQu commented Oct 11, 2021

@endolith

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.

Regards

@endolith
Copy link
Author

@ThomQu

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.

@endolith
Copy link
Author

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 * ((np.prod(2*pi*1000 - z)/
              np.prod(2*pi*1000 - p)) *
             (np.prod(1 - exp(p*T - 2*pi*1000*T)))/
              np.prod(1 - 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 https://www.khabdha.org/wp-content/uploads/2008/03/optimizing-the-magnitude-response-of-mzt-filters-mzti-2007.pdf

@endolith
Copy link
Author

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?

@endolith
Copy link
Author

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