Instantly share code, notes, and snippets.

# endolith/A_weighting.py

Last active May 28, 2024 16:55
Show Gist options
• Save endolith/148112 to your computer and use it in GitHub Desktop.
A-weighting audio files in Python
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #!/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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 commented Oct 14, 2022

so maybe I should just make an installable weighting filter library that does nothing else

and maybe include bilinear as an option so people can compare with the results of other software?

### endolith commented Oct 14, 2022

FYI, the MZTi version is accurate down to about 2.5 kHz:

### rickou commented Aug 9, 2023

@endolith your graph are very interresting !
how do you graph the type 0 limits ?

### rickou commented Aug 10, 2023

Thank you endolith !

i just found the table from standard, and i'm little bit surprised because the values are not exactly the same.

do you still have the original table ? (i do not have the standard at my side, so i use google to find parts of it 🙃)
thank a lot !

### rickou commented Aug 10, 2023

@endolith one more question about MZTi, have you push this improvement somewhere ? (i haven't found it ...)
My goal is to use the resulting coefficient in a microcontroller based hardware with IIR/FIR filtering functions, do you know if MZTi coefficient usable in those cases ? (i'm beginner in signal processing)
thank you !

### endolith commented Aug 10, 2023

@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 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 commented Sep 14, 2023

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