9
$\begingroup$

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

enter image description here

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.

Improve this question
$\endgroup$
9
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Dec 6, 2016 at 0:12
  • 1
    $\begingroup$ @robertbristow-johnson ah! nice explanation! thanks! $\endgroup$ Commented Dec 6, 2016 at 9:26
  • 1
    $\begingroup$ Another simple solution to accurate measurements is to just upsample the data first and then apply the bilinearized filter. $\endgroup$ Commented May 10, 2017 at 15:45
  • 2
    $\begingroup$ I have an implementation using MZTi ("matched z transform, improved") that works pretty well, for arbitrary fs. I need to clean it up and post it. Remind me if I don't $\endgroup$ Commented Jul 11, 2019 at 17:56
  • 1
    $\begingroup$ @user2561747 Why is this suddenly popular in the last few days? You all in the same class or something? :) gist.github.com/endolith/148112#gistcomment-3920822 $\endgroup$ Commented Oct 15, 2021 at 15:34

7 Answers 7

Reset to default
6
$\begingroup$

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.

enter image description here

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;
Improve this answer
$\endgroup$
3
  • 3
    $\begingroup$ 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? $\endgroup$ Commented Dec 8, 2016 at 13:23
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Dec 8, 2016 at 18:01
  • 1
    $\begingroup$ @MattL. Would be very grateful if you could share your design method. $\endgroup$ Commented Oct 10, 2021 at 22:55
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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? $\endgroup$ Commented Jul 10, 2019 at 21:07
2
$\begingroup$

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:

enter image description here

(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): enter image description here

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Could you elaborate on the MIM method used for the LPF? And possibly provide a small snippet of sample code? $\endgroup$ Commented Feb 24, 2020 at 16:28
  • 1
    $\begingroup$ Neat approach - do you know how well this method fares at lower sample rates? $\endgroup$ Commented Feb 24, 2020 at 16:28
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Feb 24, 2020 at 19:09
1
$\begingroup$

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};
Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Any chance of a plot of the actual amplitude response? $\endgroup$ Commented Dec 6, 2016 at 13:27
0
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Beware: this FDLS implementation appears to incorrectly call scipy freqz in places where it should call freqs. $\endgroup$ Commented Feb 25, 2019 at 21:50
0
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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()
Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ 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). $\endgroup$ Commented Jul 10, 2019 at 20:31
  • $\begingroup$ 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/… $\endgroup$ Commented Jul 10, 2019 at 20:31
  • $\begingroup$ 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 $\endgroup$ Commented Jul 10, 2019 at 20:44
  • $\begingroup$ Thanks for commenting! Indeed I forgot to append the sample to the test signal. Fixed that. $\endgroup$ Commented Jul 11, 2019 at 8:35

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.