Create a gist now

Instantly share code, notes, and snippets.

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

This comment has been minimized.

Show comment
Hide comment
@Zeokat

Zeokat Mar 5, 2014

Zeokat thanks for the code, will test it.

Zeokat commented Mar 5, 2014

Zeokat thanks for the code, will test it.

@mankind-nil

This comment has been minimized.

Show comment
Hide comment
@mankind-nil

mankind-nil Jan 13, 2015

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

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

@acuti

This comment has been minimized.

Show comment
Hide comment
@acuti

acuti 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.

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

This comment has been minimized.

Show comment
Hide comment
@frgos2

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

`

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)

`

@endolith

This comment has been minimized.

Show comment
Hide comment
Owner

endolith commented Mar 22, 2018

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