A-Weighting filter

[juha_p]

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

As separate filters:

LPF (1st order)

Code: Select all

0.495742728952490   0.201393367086572
1.000000000000000  -0.304426466162266

LPF (2nd order)

Code: Select all

0.495742728952490   0.602608294382985   0.158790522194091
1.000000000000000   0.504894384885179  -0.252335432834296

HPF (4th order, same for both)

Code: Select all

1.18767698102009  -4.75070792408037   7.12606188612055  -4.75070792408037   1.18767698102009
1.000000000000000  -3.888553450107600   5.667527095167571  -3.669386132781579   0.890412497063071

As combined filters:
5th order:

Code: Select all

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

6th order:

Code: Select all

[[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]

Both gives quite good response:

https://i.postimg.cc/BZN4gRNf/e0Ow9.png

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

I noticed that compared to analog model, "peak" around 2.5kHz is little off but, supposing, it is possible to fix by changing the dc gains of HPF and LPF a bit. There also seems to be something going on at below 2 Hz (maybe because of using c2d(.. 'matched')) but, as the magnitude level is below -120dB it shouldn't be an issue ... (maybe mapping using BLT could improve that ...?).

As I'm noobie with DSP/Matlab/Octave, how can one find error (curve) against analog model ?

EDIT: Responses for low sample rates (4, 8, 16 and 32 kHz):

https://i.postimg.cc/VLfL7xnM/9tNJa.png

[hugoderwolf]

May I ask why you are so extremely obsessive over accurately matching these weighting curves? I mean they are kind of an arbitrary historic snapshot of what someone at some point in time thought was a good idea for "measuring" (rather estimating) perceived loudness (and implementable with available parts and technology). It's not like they contain any kind of ultimate truth of measurement. In fact, they're far from it. I'd really like to understand the value of this.

[mystran]

There also seems to be something going on at below 2 Hz (maybe because of using c2d(.. 'matched')) but, as the magnitude level is below -120dB it shouldn't be an issue ... (maybe mapping using BLT could improve that ...?).

Looks like a numerical issue. It could be a measurement error (as suggested by the noise, say if the measurement is done from a truncated impulse response) or it could be an artifact of the numerical problems with a direct-form filter (which you seriously don't want to use beyond 2nd order if you can help it; unfortunately by the time you have the direct form coefficients, some of the damage is already done).

[mystran]

May I ask why you are so extremely obsessive over accurately matching these weighting curves? I mean they are kind of an arbitrary historic snapshot of what someone at some point in time thought was a good idea for "measuring" (rather estimating) perceived loudness (and implementable with available parts and technology). It's not like they contain any kind of ultimate truth of measurement. In fact, they're far from it. I'd really like to understand the value of this.

They are somewhat reasonable estimates for measuring noise-levels and you might want to match them somewhat closely if you want your measurements to match those of another measurement tool. That said, the standards define tolerances for these curves so a bit of variation is to be expected [but don't ask me about the specifics of said tolerances, I don't care enough to look it up].

That said, I don't think these curves have much practical value in software. The main problem with regards to music in particular is that the loudness sensitivity is apparently somewhat different depending on whether you are looking at noise or pure tones, with music typically coming closer to the latter. For measuring musical signals, I would rather suggest either no weighting, straight 3dB/octave weighting, or ITU-R BS.1770 / EBU R128 loudness weighting (which is also problematic, but at least it's a standard designed for measuring music), depending on exactly what you are trying to do.

[juha_p]

May I ask why you are so extremely obsessive over accurately matching these weighting curves? ...
I'd really like to understand the value of this.

Just because of I like to do so ...

[juha_p]

Looks like a numerical issue. It could be a measurement error (as suggested by the noise, say if the measurement is done from a truncated impulse response) or it could be an artifact of the numerical problems with a direct-form filter (which you seriously don't want to use beyond 2nd order if you can help it; unfortunately by the time you have the direct form coefficients, some of the damage is already done).

Thanks,
I have to dig this more next weekend.

[juha_p]

OK, digged a bit and found out that using 1st order s-domain models for all HP filters + using 'tustin' method on c2d() call improves the LF response but, still, higher the sample rate ... worse the situation.

Q: How can gain for a combined z-domain filter (FLT=LP*HP1*HP1*HP2*HP3) be matched to be 0dB at 1000Hz (Octave/MATLAB)? I did found one solution but, as it is done already at s-domain stage so, it's not usable in this implementation... .

[mystran]

Tustin's method is just another name for BLT.

edit: Also if higher sampling rate makes the low-frequency response worse, this typically suggests numerical issues, especially with direct form filters.

[juha_p]

Tustin's method is just another name for BLT.

Of course.

[mystran]

Q: How can gain for a combined z-domain filter (FLT=LP*HP1*HP1*HP2*HP3) be matched to be 0dB at 1000Hz (Octave/MATLAB)? I did found one solution but, as it is done already at s-domain stage so, it's not usable in this implementation... .

You can simply compute the response (complex or amplitude) of each of the filters at the particular frequency, then multiply them together. This will get you the combined response, which you can then divide out (ie. set the "global" gain to the reciprocal). Alternatively just normalise each of the parts separately.

[juha_p]

You can simply compute the response (complex or amplitude) of each of the filters at the particular frequency, then multiply them together. This will get you the combined response, which you can then divide out (ie. set the "global" gain to the reciprocal). Alternatively just normalise each of the parts separately.

If I normalize each filter at s-domain level, that results ~-1.98 dB off.

[mystran]

You can simply compute the response (complex or amplitude) of each of the filters at the particular frequency, then multiply them together. This will get you the combined response, which you can then divide out (ie. set the "global" gain to the reciprocal). Alternatively just normalise each of the parts separately.

If I normalize each filter at s-domain level, that results ~-1.98 dB off.

If you want to normalize the final response, then you should normalize the final filters after any transform/fitting process. Please keep in mind though that in some cases this can increase the overall error.

[juha_p]

OK, here's the final Octave code to produce 5th order IIR filter coefficients for a A-Weighting filter:

Code: Select all

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

As an alternative for building HPF filters from analog model, one can use this function:

Code: Select all

% This 1st order (Butterworth type) HPF coefficient calculation function 
% uses low order approximations and is meant to be used for cut off frequencies 
% below 1 kHz
% Approximation error is less than ±5e-05 in (frequency) range [1 Hz; 1000 Hz]
%
% Usage: [b,a] = HPF1_approx(44100, 350, 2)
%
% Brought to you by jiiteepee@yahoo.se (02/2020)
% ------------------------------------------------------------------------------
function [b, a] = HPF1_approx(fs, fc, order)
% fs    = sample rate in Hz
% fc    = cut off frequency in Hz
% order = approximation order; 3 = 3rd order, otherwise 2nd order

x = 2.0 * pi * fc/fs;

switch order
  case 3 % 3rd order
    sx = -5.382643908e-8 + x * (1.00000755295 + x * (-1.627621479421707e-4 + ...
          x * (-0.1657084756226908)));
    cx = 1.99999946747854 + x * (8.3526171644607921e-5 + x * (-0.50208137392503 + ...
          x * 1.664302325330385905e-2));
  otherwise % 2nd order
    sx = -4.205541129472213e-5 + x * (1.003746378622 + x * (-4.991453468087422e-2));
    cx = 2.00000427449401 + x * (-3.1321682298812e-4 + x * (-0.4969816314888122));
end

% return
a = [cx+sx sx-cx];
b = [cx -cx];

Which is used as:

Code: Select all

[b, a] = HPF1_approx(Fs,f1, 3);
HP1=tf(b, a, 1/Fs);
...