DSPRelated.com
Forums

A-Weighting filter

Started by jtp_1960 4 years ago3 replieslatest reply 4 years ago311 views
A-Weighting filter implemented here

uses Bilinear transformation method, which has some issues at high frequency area. 

Here's one alternative I made using Octave:

<code>% 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: <a href="https://www.dsprelated.com/showcode/214.php">https://www.dsprelated.com/showcode/214.php</a>
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;
</code>

which produces quite accurate magnitude response already for 5th order filter implementation:

aweighting2_56538.png

and works well for low samplerates as well (Fs = 4, 8, 16 and 32 kHz):


aweightlowfs_1135.png



[ - ]
Reply by aclarkFebruary 29, 2020

I gave a paper that addressed this problem at the 2010 comp.dsp conference.

I have it in pdf format, but I am not sure how to submit it.

Here are some of the ideas.


If you split the high pass portion of the A weight filter from the low pass part, then you can use BZT for the high pass section and get a very good fit for this part. This leaves you with the double real poles at 12.2K. This portion is fitted with a magnitude squared approach.

Therefore rewrite

HA(s) = H123(s) H4(s)

HA(s) = [s4/[(s+w1)2(s+w2)(s+w3)] *

w42 /(s+w4)2

H123(s) -> HBZT123(z) Use BZT

H4(s) = [w4 /(s+w4)] * [w4 /(s+w4)]

Use two second order HFIR (z) where

1st section uses w0 = 0, w1 = wp, w2 = p-wp

2nd section uses w0 = 0, w1 = p/3, w2 = 2p/3

* Using two different matching criteria spreads the error

* w0, w1, w2 are not the A weight definitions

* wp = 76618.5439, the pole location Magnitude squared


If HFIR(z) = b0 + b1z-1 + b2z-2 then

|HFIR(z)|2= b02 + b12 + b22 + 2b0b1cos w +

2b0b2cos 2w + 2b1b2cos w

We now compute the filter by matching three points. This is difficult to factor for arbitrary w, but the above choices are easier to solve.



Part of this idea comes from this paper:

David Gunness & Ojas Chauhan, Optimizing the magnitude response of matched Z transform filters (MZTi) for Loudspeaker Equalization, Proceedings AES 32nd International Conference 21-23 Sept 2007


Al Clark

[ - ]
Reply by jtp_1960February 29, 2020

That's much easier to implement compared to Magnitude Invariance Method!

What comes to submitting PDF ... doesn't that file inserting (paperclip symbol) feature support PDF?

[ - ]
Reply by aclarkFebruary 29, 2020