# Design of a digital A-weighting filter with arbitrary sample rate

I want to A-weight a time series with arbitrary sample rate.

An analog A-weighting filter is defined exactly by IEC 61672-1. But there's no definition for a digital filter. One method is to use the bilinear transform (BLT) to convert the analog filter to the digital filter (as done here Applying A-weighting). However this method suffers from extreme warping near nyquist (even when the analog poles/zeros are pre-warped):

Figure 1: A-weighting frequency response comparison where the sample rate is $25600\textrm{ Hz}$.

Instead I'm thinking of using an algorithm than can design a digital IIR filter with arbitrary frequency response and plugging in the frequency response of the analog A-weighting filter.

• Is this a good approach?
• If so, is there a particular algorithm that would be well suited for this?

I've looked into MATLAB's yulewalk but I would need a corresponding Python implementation to try out. I've also come across Berchin's FDLS method in a few places, like this question for instance, but all of the links appear to be broken.

• personally, I don't really see the problem with your digital filters – yes, there'll be phase warping, but if you look at the analog filter, it leads to aliasing – so your filter damages the signal "less". – Marcus Müller Dec 5 '16 at 23:47
• the problem @MarcusMüller is that, when using the bilinear transform, what the analog filter does at $\infty$ is what the digital filter does at $\pi$ or Nyquist. every frequency specification of the original filter can be pre-warped (so that it comes out correctly after applying bilinear transform) but the behavior at Nyquist cannot be pre-warped. so then you have to do the Orfanidis thing and design the analog filter to do at $\infty$ what you want the digital filter to do at Nyquist. – robert bristow-johnson Dec 6 '16 at 0:12
• @robertbristow-johnson ah! nice explanation! thanks! – Marcus Müller Dec 6 '16 at 9:26
• I'm not sure whether it's the best fit for this application or not, but here is another question on Berchin's FDLS that might be enough for you to implement it and try. – Jason R Dec 8 '16 at 18:12
• Another simple solution to accurate measurements is to just upsample the data first and then apply the bilinearized filter. – endolith May 10 '17 at 15:45

It's a common misconception that the approximation of an analog filter by a digital filter must be bad close to Nyquist. This idea might come from the ubiquity of the bilinear transform, for which this is usually indeed the case. Of course, there are certain constraints on the frequency response of discrete-time filters at Nyquist, but they do not necessarily need to result in a bad approximation of an analog filter in that frequency range. The quality of the approximation close to Nyquist depends on several factors, among which are the properties of the frequency response of the analog filter, and the fact whether only the magnitude or also the phase of the analog filter need to be approximated.

I've designed a 6th order IIR filter approximating the frequency response of an analog A-weighting filter, as defined here:

$$H(s)=\frac{k\cdot s^4}{(s+129.4)^2(s+676.7)(s+4636)(s+76655)^2}\tag{1}$$

with $k=7.39705×10^9$.

I chose a sampling frequency of $48$ kHz. The design procedure is a heuristic iterative procedure I came up with some time ago. It's a least squares approximation based on the equation error method, and I might write up all the details some day.

Below is a plot of the design result. Note that both plots go up to Nyquist ($24$ kHz). The left figure shows that one can't see any difference between the logarithmic plots of the analog and the digital frequency responses. The right-hand figure shows the approximation error, defined as the absolute value of the difference of the magnitude responses. The error shows a typical least squares behavior.

You can check out the filter yourself. Here are the coefficients:

b =

0.169994948147430
0.280415310498794
-1.120574766348363
0.131562559965936
0.974153561246036
-0.282740857326553
-0.152810756202003

a =

1.00000000000000000
-2.12979364760736134
0.42996125885751674
1.62132698199721426
-0.96669962900852902
0.00121015844426781
0.04400300696788968

• Nice answer. If you could give slightly more detail on how you generated the coefficients, then this would be an even better answer, since the OP was asking for designing the filter for an arbitrary sample rate. Is your technique similar to Berchin's frequency-domain least-squares (FDLS) method? – Jason R Dec 8 '16 at 13:23
• @JasonR: Thanks for your comment, I will add some more information as soon as I have the time. As for FDLS, I don't know all the details, and I haven't yet seen an authorized implementation, but I know that FDLS needs a complex desired response (i.e., magnitude and phase), whereas here we're only interested in the magnitude. – Matt L. Dec 8 '16 at 18:01

Depends on how fancy you need to get. A good implementation at 48 kHz is the following

b =    [0.234301792299513  -0.468603584599026  -0.234301792299513  ...
0.937207169198054  -0.234301792299515  -0.468603584599025   0.234301792299513];
a =    [1.000000000000000  -4.113043408775871   6.553121752655047 ...
-4.990849294163381   1.785737302937573  -0.246190595319487   0.011224250033231];


If you need it at a different sample rate, you can simply calculate the impulse response of the IIR filter at 48 kHz, truncate and resample to whatever sample rate you want. At 48 kHz, 256 taps gets you down to about 50 Hz. 1024 taps down to 20 Hz (if you care) and 2048 is as good as you would ever need it.

Resampling will always impact the frequency response at the Nyquist frequency but that's a fact of life for all signal processing task.

Original filter coefficients courtesy of Dr Ir Christophe Couvreur Faculte Polytechnique de Mons
Details removed for spam protection, let me know if you want contact info.

• When I initially read your reply (years ago), I mistakenly thought you were suggesting to resample the input signal. Now (with more experience) I understand you're suggesting to resample the impulse response of an FIR filter (much shorter than the input signal). This technique sounds very promising; can you provide a little more details, such as which window you used for truncation and which resampling method? – user2561747 Jul 10 '19 at 21:07

Here are coefficients for 5th and 6th order filters (fs=48kHz) I prepared using separate HPF (c2d)) and LPF (MIM):

4th order HPF, 1st order LPF:

[0.588402730019084 -2.114578549207340 2.574286896638522 -0.919416694862366 -0.367726753456895 0.239032370868995]
[1.0 -4.193437345479311 6.853084418397484 -5.397323870680988 2.009149234848459 -0.271472430592242]


4th order HPF, 2nd order LPF:

[[0.5884027300190836 -1.6384145461593100 0.8580879826180272 1.1837389307393180 -1.1416401766192728 -0.0386320187694699 0.1884570981716236]
[1.0 -3.384188885614779 3.453610828634393 0.171625555633360 -2.392296872830998 1.376227862809258 -0.224978476938553]


and a plot showing comparison against their analog model:

(green = analog model, red = 5th order filter, black = 6th order filter, small window shows the difference of responses when approaching the Nyqvist frequency)

EDIT2:

Here's Octave source code for to build the filter:

% Octave packages -------------------------------

% other requirements:
% Magnitude Invariance method (MIM) -implementation ( one implementation can be found from https://soar.wichita.edu/handle/10057/1564 )

clf;
format long;

%Sampling Rate
Fs = 44100;

% A-weighting filter frequencies according to IEC/CD 1672.
% Source: https://www.dsprelated.com/showcode/214.php
f1 = 20.598997;
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;

% Analog model -----------------------------------------------------
A1000 = 1.9997;
NUM = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
DEN = conv(conv(DEN,[1 2*pi*f3]),[1 2*pi*f2]);

Ds = tf(NUM, DEN);
[Ab, Aa, T] = tfdata(Ds, 'v');

% Digital filter  --------------------------------------------------
%LPF1 (MIM method)
lporder = 1; % 1 = 5th order A-Weighting filter, 2= 6th order ....
w0 = 2 * pi * f4;
bC = 1;
aC = [1 w0];
AD0 = tf(bC, aC);
LP2 = c2dn(AD0^2, 1/Fs, 'mim', lporder, 4096^2, 'lowpass');

% HPF (BLT method)
% HPF1
w0 = 2 * pi * f1;
bC = [w0 0];
aC = [1 w0];
AD1 = tf(bC, aC);
HP1 = c2d(AD1^2, 1/Fs, 'tustin');

% HPF2
w0 = 2 * pi * f2;
bC = [w0 0];
aC = [1 w0];
AD2 = tf(bC, aC);
HP2 = c2d(AD2, 1/Fs, 'tustin');

% HPF3
w0 = 2 * pi * f3;
bC = [w0 0];
aC = [1 w0];
AD3 = tf(bC, aC);
HP3 = c2d(AD3, 1/Fs, 'tustin');

% Combine filters
FLT = (LP2*HP1*HP2*HP3);

GAINoffset = abs(freqresp(FLT,1000*2*pi));
FLT = FLT/GAINoffset;

% get coefficients
[ad, bd, T] = tfdata(FLT, 'v');
bd

% plot
fs2 = Fs/2;
nf = logspace(0, 5, fs2);

[mag, pha] = bode(Ds,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'g', 'linewidth', 4.0, 'linestyle', '-');

hold on;

[mag, pha] = bode(FLT,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'k', 'linewidth', 2.0, 'linestyle', '--');
title('A-Weighting filter');
xlabel('Hz');ylabel('dB');
axis([1 fs2+5000 -200 10]);

legend('Analog model', 'Digital (MIM+BLT)', 'location', 'southeast');
grid on;


I have not measured the error curve (dunno how to calculate it against analog model).

EDIT1: Plot showing HF response for low sample rates (4, 8, 16, 32 kHz):

• Could you elaborate on the MIM method used for the LPF? And possibly provide a small snippet of sample code? – user2561747 Feb 24 at 16:28
• Neat approach - do you know how well this method fares at lower sample rates? – user2561747 Feb 24 at 16:28
• Works well @ lower sample rates ... actually that low frequency issue gets worse at higher sample rates ... . Regarding MIM method, paper (thesis) with matlab code can be found from here: soar.wichita.edu/handle/10057/1564 . Should work for LPF without modifications. – Juha P Feb 24 at 19:09

I made my own matlab routine a long time ago to give me IIR filter coefficients for A-weighting. From a now-dead link, I got the transfer function of the form:

                              ka*s^4
------------------------------------------------
(s+129.4)^2 * (s+676.7) * (s+4636) * (s+76655)^2


From this starting point, I used matlab to build a discrete version of the filter. It's ugly, but here it is:

 function [b,a]=makeAweightingFilter(fs);

%build the continuous-time transfer function
num=poly([0 0 0 0]);    %numerator of transfer function (TF)
den=poly(-129.4336);   %one of the poles of the TF
den=conv(den,poly(-129.4336));   %another pole of the TF
den=conv(den,poly(-676.6991));   %another pole of the TF
den=conv(den,poly(-4636.3624));  %another pole of the TF
den=conv(den,poly(-76654.86075));  %another pole of the TF
den=conv(den,poly(-76654.86075));  %another pole of the TF
a_weight=tf(num,den);  %form the transfer function

%Scale transfer function to yield 0dB response at 1000 Hz
offset = abs(freqresp(a_weight,1000*2*pi));
num=1/offset*poly([0 0 0 0]); %scale the numerator
a_weight=tf(num,den); %re-form the transfer function

%Convert to discrete time system
a_w_discrete = c2d(a_weight,1/fs,'zoh');

%get filter coefficients from the numerator
%and denominator of the discrete transfer function
b=a_w_discrete.num{1};
a=a_w_discrete.den{1};

• Any chance of a plot of the actual amplitude response? – jojek Dec 6 '16 at 13:27

There is a Python implementation of Berchin FDLS here: https://github.com/happycube/ld-decode/blob/a9418b96a492ca6b6d6b0bab33a244b371e5e348/fdls.py

This could be used to yield an IIR which approximates the transfer function of the analog filters.

• Beware: this FDLS implementation appears to incorrectly call scipy freqz in places where it should call freqs. – user2561747 Feb 25 '19 at 21:50

There is a nice one-shot pseudo-least squares method described in "A Generalization of the Biquadratic Parametric Equalizer" by Knud Bank Christensen that I have succesfully used for a few things, including designing an A-weighting filter.

http://www.aes.org/e-lib/browse.cfm?elib=12429

If you have Matlab, the actual filter design bit is implemented in the invfreqz() function, then you do some pre and post frequency warping. From memory though I think the mothod is described well enough in the appendix to implement in Python without much trouble.

I had to play with parameters a bit as I was implementing the final filter in fixed point and it is easy to get very large coefficients out of this. But with a bit of juggling I got something that was good enough for me (but admittendly not as accurate as Matt Ls solution). If you're working in floating point you won't have to worry about this though.

Hinted by the answer from @chipaudette I made a Python implementation using the s-domain transfer function (https://en.wikipedia.org/wiki/A-weighting).

# python 3.7
import scipy.signal as signal
import math as math
import matplotlib.pyplot as plt
import array as arr

# for discretization
sr = 48000.0
# enter the zeros, poles and scaling as in s-domain transfer function from Wikipedia
z,p,k = signal.bilinear_zpk([0.0, 0.0, 0.0, 0.0], [-129.4, -129.4, -676.7, -4636.0, -76655.0, -76655.0], 7.39705E9, sr)
# get the filter coefficients
b,a = signal.zpk2tf(z,p,k) #delivers 'ba' output
# generate test signal of 1 second
test_signal = arr.array('d',[])
for x in range(int(sr)):
v = math.sin(x * 1000 * 2.0 * math.pi / sr) #1000hz test signal
test_signal.append(v)
# apply filterapply our a-filter
filtered = signal.lfilter(b, a, test_signal)


If you want to see the frequency response of the filter, that can be done easily in the continuous domain.

# python 3.7
import scipy.signal as signal
import math as math
import matplotlib.pyplot as plt
import array as arr
# make the system continuous
system = signal.ZerosPolesGain([0.0, 0.0, 0.0, 0.0], [-129.4, -129.4, -676.7, -4636.0, -76655.0, -76655.0], 7.39705E9)
# plot the frequency response
w, mag, phase = signal.bode(system)
plt.figure()
plt.semilogx(w, mag)
plt.show()

• Hi, thanks for posting fully reproducible code! This is the bilinear transform method mentioned in the question, which suffers from frequency warping near nyquist (which can be a problem at lower sample rates). – user2561747 Jul 10 '19 at 20:31
• Also be warned: it's not meaningful to look at the frequency response of the continuous filter, which does not match that of the digital filter (see the plot in the question - "Analog" is the response of the continuous filter, and "BLT No Prewarp" is a digital filter constructed from the continuous filter using the bilinear transform as you've done). When using digital filters it's important to look at the response of the actual digital filter being used, which can be done with docs.scipy.org/doc/scipy/reference/generated/… – user2561747 Jul 10 '19 at 20:31
• The test signal code doesn't work unfortunately, since v is getting overwritten before being appended to anything. I highly, highly recommend using numpy.org for working with arrays, which you should already have installed as it's a dependency of scipy. It can simplify and speed up generation (try generating an hour long signal and compare performance ;) of a test signal: import numpy as np sr = 48000.0 # Hz test_signal = np.sin(np.arange(sr) * 1000 * 2 * np.pi / sr) # 1000 Hz test signal – user2561747 Jul 10 '19 at 20:44
• Thanks for commenting! Indeed I forgot to append the sample to the test signal. Fixed that. – Lodewijk Loos Jul 11 '19 at 8:35