A-weighting audio files in Python

readme.md

You should use this version instead; it's better.

Note that this uses a bilinear transform and so is not accurate at high frequencies.

Apply an A-weighting filter to a sound stored as a NumPy array.

Use Audiolab or other module to import .wav or .flac files, for example. http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/audiolab/

Translated from MATLAB script (BSD license) at: http://www.mathworks.com/matlabcentral/fileexchange/69

A_weighting.py

#!/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)

analyzer.py

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

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.

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)