C-Weighting filter

[juha_p]

Trying to implement "accurate" C-weighting filterby using other than bi-linear transformation method.

What I've achieved so far is 10th order filter by using Magnitude Invariance Method (MIM (Matlab implementation)):

https://i.postimg.cc/rsw8h6X7/CcNyE.png

but, as the method has an issue with HPF --> this is the lowest order implementation I can accept the result of.

I know it does not need to be that high order filter to get accurate response but, so far I have not found a good solution.

Accurate HPF response can be done with low order BLT but the LPF part needs tweaking (low order MIM can do this well enough). So, is it possible to combine a BLT type filter with a non-BLT filter (using Matlab/Octave) or, ... is it the only way to go by trying to improve BLT type filter response ?

[matt42]

Why not some kind of MZT with curve fitted zeroes? Your curve looks fairly simple probably you would get extremely close with a first order low pass and a first order high pass. 10th order seems like overkill for that curve

[juha_p]

Aren't analog c-weighting filter done by using two 1st order HP filters and two 1st order LP filters connected in series ... ?

I was thinking if combing z-domain filters is possible ... this way I could make HP filtering by combining two 1st order as biquad (BLT) implementation and LP filtering using two 1st order LP filters as 3rd order (MIM) implementation ... and combine these to one filter but is it possible to do?

[matt42]

Aren't analog c-weighting filter done by using two 1st order HP filters and two 1st order LP filters connected in series ... ?

Could well be I don't know what the transfer function is. But my point is that the curve is very simple and would be easily matched with a modified MZT filter.

If you want to combine BLT LPF with MIM high pass I don't see what would stop you doing that. I would just run them as seperate filters. Otherwise if you want a combined structure you'd just multiply the z-domain transfer functions.

[mystran]

I was thinking if combing z-domain filters is possible ... this way I could make HP filtering by combining two 1st order as biquad (BLT) implementation and LP filtering using two 1st order LP filters as 3rd order (MIM) implementation ... and combine these to one filter but is it possible to do?

You would usually want to factor higher order filters into biquad sections anyway (for numerical stability), so why not just put the two filters in series?

[juha_p]

Yes, using 2nd order sections is more secure but is that OK with metering use?

Actually, I was thinking, when you want better accuracy, it has to be a higher order implementation than what the original is.

Example:
44100 Hz implementation using Matlab example behind the link in my opening post:

s-domain coefficients:
[5.9124e+09 1.0000e-32 1.0000e-32]
[1.0000e+00 1.5350e+05 5.9101e+09 1.5221e+12 9.8338e+13]

BLT
[0.21701 0.00000 -0.43402 0.00000 0.21701]
[1.0000000 -2.1346750 1.2793335 -0.1495598 0.0049087]

My implementation:

MIM (LPF)
s-domain coefficients:
[1]
[1.0000e+00 1.5324e+05 5.8704e+09]

[0.53249 0.20838]
[1.00000 -0.25830]

or 2nd order version
[0.53249 0.63248 0.15785]
[1.00000 0.53815 -0.21811]


BLT (HPF)
s-domain coefficients:
[1.2943e+02]
[1.0000e+00 2.5885e+02 1.6751e+04]

[0.99707 -1.99414 0.99707]
[1.00000 -1.99414 0.99415]

So, could a 3rd order filter be possible by just combining that 1st order LPF with HPF?

[mda]

I got a good result with 12 dB highpass, 12 dB lowpass, and two Nyquist resonators to fix the HF roll-off. So effectively 6-pole.

For the resonators I mean a one-pole filter

Code: Select all

 out += a * (in - out)

where a > 1.0

[mystran]

Yes, using 2nd order sections is more secure but is that OK with metering use?

If we ignore numerical precision (for just a moment), then any linear time-invariant filter is completely defined by it's transfer function poles and zeroes and how you choose to implement said filter doesn't have any effect on the final result. As long as the poles and zeroes are the same, the result is the same and everything else is an implementation detail.

In practice, when we work in finite-precision, we typically want to try and minimise any numerical errors, because the actual response of a poorly behaved implementation might be quite far from the intended theoretical response. Now 1st order sections would be the most stable, but we need at least 2nd order sections to realise complex-conjugate pairs of poles or zeroes in a real-valued filter and their precision is still quite fine (where as especially for direct form filters, it gets rapidly worse as the order increases), so it usually makes sense to just cascade 2nd order sections (with an extra 1st order section for odd-order filters).

The important point to understand though is that in theory this has no effect on the result, while in practice it just gives us a result closer to the theoretical ideal we were trying to get.

ps. The implementation structure does become important when you have a filter that is either time-varying and/or non-linear, but this is almost certainly not what you want for measurement purposes.

pps. I'd also like to note that if you have two 1st-order filters, then multiplying the transfer functions together to get a single 2nd-order filter is generally fine... but you probably should avoid going any higher order like this, because the numerical precision really does become a problem much faster than you'd think. In fact, it even typically makes sense to go the extra mile and write your filter design functions in such a way that they can work with biquad sections directly, because trying to factor a numerical higher-order transfer function back into 2nd order sections also tends to be quite numerically unstable.

[juha_p]

Ok,

I think implementation using two biquad filters (i.e. 4th order system) gives quite good result already (thought the 1st order LPF would be usable as well).

Here's the plot showing filters done using data listed in my previous post.

https://i.postimg.cc/0NzYqGFC/cweighting-2xbiquads.png

[juha_p]

...
Your curve looks fairly simple probably you would get extremely close with a first order low pass and a first order high pass.
...
I don't know what the transfer function is. But my point is that the curve is very simple and would be easily matched with a modified MZT filter.
...

MIM method looks like kind of MZT/IIM variant but it does not work well with HP.?
It seems that MIM does very good job in discretization the original s-domain (LPF part) filter as a 1st order filter.

[mystran]

...
Your curve looks fairly simple probably you would get extremely close with a first order low pass and a first order high pass.
...
I don't know what the transfer function is. But my point is that the curve is very simple and would be easily matched with a modified MZT filter.
...

MIM method looks like kind of MZT/IIM variant but it does not work well with HP.?
It seems that MIM does very good job in discretization the original s-domain (LPF part) filter as a 1st order filter.

No personal experience with MIM (though I did look at a very short paper about it), but my experience with other fitting methods (eg. MZTi and pole-zero mapping with a separate corrector FIR) suggests that you can often make things work, by mapping the problematic s-domain zeroes directly, then essentially adding their response backwards to the error term of the remaining filter. Another trick that could work (given that I guess the main issue here is the cepstrum computation) is to move the zeroes very slightly from the imaginary axis to allow for the logarithm to stay finite (eg. make the DC response -200dB or something, rather than a canonical zero). In any case, the general theme here is that you can usually make things work, but you might have to adjust the whole fitting method slightly.

That said, BLT is generally perfectly fine if the poles and zeroes of the system are at a low frequency (eg. below fs/6 or something is a good rule of thumb, though it depends somewhat on the response), which is certainly the case for the HP part of your target response, so I personally would just use the result you plotted and be done with it.

[juha_p]

EDIDED:

OK, here's a here's the final Octave code to produce 3rd order IIR filter coefficients for a C-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 )


format long;
clf;

%Sampling Rate
Fs = 44100;

%Analog C-weighting filter according to IEC/CD 1672.
% Source: https://www.dsprelated.com/showcode/215.php
f1 = 20.598997; 
f4 = 12194.217;

C1000 = 0.0619; %% adjust gain 

NUM = [ (2*pi*f4)^2*(10^(C1000/20)) 1e-32 1e-32];
DEN_4 = [1 +4*pi*f4 (2*pi*f4)^2];
DEN_1 = [1 +4*pi*f1 (2*pi*f1)^2];
DEN = conv(DEN_4,DEN_1);

Ds = tf(NUM, DEN);
[b, a, T] = tfdata(Ds, 'v');

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

% HPF1
w0 = 2 * pi * f1;
bC = [w0 err];
aC = [1 w0];
AD = tf(bC, aC);
HP2 = c2d(AD^2, 1/Fs, 'matched');

FLT = LP2*HP2; 
offset = abs(freqresp(FLT,1000*2*pi));
FLT = FLT/offset;

[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', 5.0, 'linestyle', '-');
hold on;

[mag, pha] = bode(FLT,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'k', 'linewidth', 2.0, 'linestyle', '--');
        
title('C-Weighting Filter');
legend('Analog model', 'Digital (MIM+BLT)', 'location', 'southeast');
xlabel('Hz');ylabel('dB');
axis([1 fs2 -55 5]);
grid on;

and couple plots:

https://i.postimg.cc/8CHHqkxH/cweighting.png

and above 5 kHz
https://i.postimg.cc/xdY5xNp8/cweign.png

EDIT:

As an alternative for building HPF filter 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);