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.

Zeokat commented Mar 5, 2014

Zeokat thanks for the code, will test it.

@mankind-nil

This comment has been minimized.

mankind-nil commented 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.

@acuti

This comment has been minimized.

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.

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.

Owner

endolith commented Mar 22, 2018

@shangyan0517

This comment has been minimized.

shangyan0517 commented Aug 2, 2018

Dear all, I am new to digital signal processing. I am wondering, in order to get A weighted sound pressure level, should I first apply this function to the wav data read by functions like "scipy.io.wavfile.read" then convert it to SPL, or first convert wav data to sound pressure level then use this function for A weighting? Thanks in advance!

@endolith

This comment has been minimized.

Owner

endolith commented Aug 2, 2018

@shangyan0517 I would recommend:

  1. Using this version instead
  2. Use PySoundFile instead of scipy.io.wavfile.read

Then you need to figure out what your dBSPL ↔ dBFS conversion factor is. Depends on sensitivity of microphone, gain of preamp, etc.

It doesn't matter if you do SPL conversion first, then A-weighting, or A-weighting, then SPL conversion. Both are LTI operations.

@shangyan0517

This comment has been minimized.

shangyan0517 commented Aug 3, 2018

@endolith Thank you so much!

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