Create a gist now

Instantly share code, notes, and snippets.

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))))

Zeokat commented Mar 5, 2014

Zeokat thanks for the code, will test it.

Thanks very much :-) The Couvreur Matlab scripts are awesome and used all over the world! Good tip about 24bit audio files too.

acuti commented Oct 25, 2017

Many thanks endolith! I am using your A_weighting.py script in some code I have started drafting, and I have included it in my repository: https://github.com/acuti/acuSLM. I hope you are ok with it. If not, please let me know.

frgos2 commented Oct 31, 2017

And here is the C-Weighting:

def C_weighting(fs):
    f1 = 20.598997 
    f4 = 12194.217
    C1000 = 0.0619

    NUMs = [(2*pi*f4)**2*(10**(C1000/20.0)),0,0]
    DENs = polymul([1,4*pi*f4,(2*pi*f4)**2.0],[1,4*pi*f1,(2*pi*f1)**2]) 

    #Use the bilinear transformation to get the digital filter. 
    return bilinear(NUMs,DENs,fs)

`

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