Skip to content

Instantly share code, notes, and snippets.

@endolith
Last active March 15, 2024 07:12
Show Gist options
  • Star 32 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save endolith/148112 to your computer and use it in GitHub Desktop.
Save endolith/148112 to your computer and use it in GitHub Desktop.
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))))
@endolith
Copy link
Author

@rickou The values I used are from https://law.resource.org/pub/us/cfr/ibr/002/ansi.s1.4.1983.pdf

I think I only posted the MZTi code above https://gist.github.com/endolith/148112?permalink_comment_id=3928477#gistcomment-3928477

and again I have qualms about it, I'm not sure if it's handling the high-frequency poles/zeros correctly.

@rickou
Copy link

rickou commented Aug 11, 2023

Thank you ! now i understand the differences..

About the MZTi, i already read your post, do you have any clues that could explain why the high frequency poles/zeros could not be good ?
What can be wrong on signal or on response ?
Have you made some test with real waveforms ?

@endolith
Copy link
Author

do you have any clues that could explain why the high frequency poles/zeros could not be good ?

MZT appeals to me for high sample rates because the A-weighting filter is defined by precise pole locations, so it feels right to transform them this way and get a similar response. The problem is that when the sampling rate is too low, they are outside the Nyquist range and it introduces error. The i stage will then correct for this as much as possible, so it's still accurate down to 2.5 kHz, but it just feels kind of wrong to me to produce the wrong response and then patch over it with this bandaid.

@rickou
Copy link

rickou commented Sep 14, 2023

hum.. i think i understand.. to be honest i'm not so comfortable with digital filtering.
i have implemented this filter in embedded software and i will made some tests in next weeks in real situation comparing my implementation and sonometer.
but my preliminary tests on bench shows nothing bad. (i'm using 20k sample rate)

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