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

Question

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.

Answer 1

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

---

Below is the Octave/Matlab code I used to design above filter. The function eqnerror.m can be found here. Note that the iteration step is purely heuristic and it might not work well for other specificiations.

fs = 48000;
nw = 500;
w = logspace( -7, pi, nw );
w = w(:);
s = 1i * w * fs;

% analog filter
% from https://en.wikipedia.org/wiki/A-weighting#Transfer_function_equivalent
ka = 7.39705e9;
Ha = ka * s.^4 ./ ( (s+129.4).^2 .* (s+676.7) .* (s+4636) .* (s+76655).^2 );

W = ones( nw, 1 );
H = Ha;
Imax = 50;
N = 6;

for k = 1:Imax
    D = abs( Ha ) .* exp( 1i * angle(H) );
    [b,a] = eqnerror( N, N, w, D, W, 30 );
    H = freqz( b, a, w );
end

% map any poles outside the unit circle into the circle
p = roots(a);
g = 1;
for k = 1 : N,
    if abs( p(k) ) > 1,
        g = g * abs( p(k) );
        p(k) = 1 / conj( p(k) );
    end
end
a = g * poly(p); a = a(:);
sc = a(1); a = a / sc; b = b / sc;

Answer 2

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.

Answer 3

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 -------------------------------
pkg load signal

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

% Adjust 0dB@1kHz
GAINoffset = abs(freqresp(FLT,1000*2*pi));
FLT = FLT/GAINoffset;

% get coefficients
[ad, bd, T] = tfdata(FLT, 'v');
ad
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):

Answer 4

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};

Answer 5

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.

Answer 6

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.

Answer 7

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