Trying to implement a digital A frequency filter

Question

I'm trying to implement a digital filter that has the frequency response shape equal to the image below:

Where i will use equation (11) to implement it with a sampling frequency of 48kHz.

The filter coefficients can be found in the same document: Where each w' follows the formula above.

So i put everything into matlab:

fs = 48000;
f1 = 20.598997;
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
w1 = 2tan(pi(f1/fs));
w2 = 2tan(pi(f2/fs));
w3 = 2tan(pi(f3/fs));
w4 = 2tan(pi(f4/fs));
%testeb2 = 2tan(pi(250/1000))2tan(pi(250/1000))(1/sqrt(2))
% Filter coefficients for the A weighting filter
a0 = 64 + (16w2w1w3) + (4w2w1w1w3) + (32w2w1w4) + (8w2w1w1w4) + (32* w1 * w3 * w4) + (16 w2  w3 * w4) + (64 w1) + (32 w2) + (32 w3) + (64 w4) + (32 w2  w1) + (8 w2  w1w1) + (16 w1w1) + (16 w2 * w1 * w3 * w4) + (4 w2  w1w1  w3 * w4) + (32 w1  w3) + (16 w2  w3) + (8 w1w1 * w3) + (64 w1  w4) + (32 w2  w4) + (32 w3  w4) + (16 w1w1 * w4) + (8 w1w1 * w3 * w4) + (16 w4w4) + (16 w4w4 * w1) + (4 w4w4 * w1w1) + (4 w4w4  w1 * w2 * w3) + (w4w4  w1w1  w2 * w3) + (8 w4w4 * w1 * w2) + (2 w4w4 * w1w1  w2) + (8 w4w4 * w2) + (8 w4w4 * w3) + (8 w4w4 * w1 * w3) + (2 w4w4 * w1w1  w3) + (4 w4w4 * w2 * w3)
a1 = -128 + (64 w2  w1 * w3) + (24* w2 * w1w1  w3) + (128 w2  w1 * w4) + (48 w2  w1w1  w4) + (128 w1  w3 * w4) + (64 w2  w3 * w4) + (64 w2  w1) + (32 w2  w1w1) + (32 w1w1) + (96 w2 * w1 * w3 * w4) + (32 w2  w1w1  w3 * w4) + (64 w1  w3) + (32 w2  w3) + (32 w1w1 * w3) + (128 w1  w4) + (64 w2  w4) + (64 w3  w4) + (64 w1w1 * w4) + (48 w1w1 * w3 * w4) + (32 w4w4) + (64 w4w4 * w1) + (24 w4w4 * w1w1) + (32 w4w4  w1 * w2 * w3) + (10 w4w4 * w1w1  w2 * w3) + (48 w4w4 * w1 * w2) + (16 w4w4 * w1w1  w2) + (32 w4w4 * w2) + (32 w4w4 * w3) + (48 w4w4 * w1 * w3) + (16 w4w4 * w1w1  w3) + (24 w4w4 * w2 * w3)
a2 = -192 + (48 w2  w1 * w3) + (52 w2  w1w1  w3) + (96 w2  w1 * w4) + (104 w2  w1w1  w4) + (96 w1  w3 * w4) + (48 w2  w3 * w4) - (320 w1) - (160 w2) - (160 w3) - (320 w4) - (96 w2  w1) + (24 w2  w1w1) - (48 w1w1) + (208 w2 * w1 * w3 * w4) + (108 w2  w1w1  w3 * w4) - (96 w1  w3) - (48 w2  w3) + (24 w1w1 * w3) - (192 w1  w4) - (96 w2  w4) - (96 w3  w4) + (48 w1w1 * w4) + (104 w1w1 * w3 * w4) - (48 w4w4) + (48 w4w4 * w1) + (52 w4w4 * w1w1) + (108 w4w4  w1 * w2 * w3) + (45 w4w4 * w1w1  w2 * w3) + (104 w4w4 * w1 * w2) + (54 w4w4 * w1w1  w2) + (24 w4w4 * w2) + (24 w4w4 * w3) + (104 w4w4 * w1 * w3) + (54 w4w4 * w1w1  w3) + (52 w4w4 * w2 * w3)
a3 = 512 - (128 w2  w1 * w3) + (32 w2  w1w1  w3) - (256 w2  w1 * w4) + (64 w2  w1w1  w4) - (256 w1  w3 * w4) - (128 w2  w3 * w4) - (256 w2  w1) - (64 w2  w1w1) - (128 w1w1) + (128 w2 * w1 * w3 * w4) + (192* w2 * w1w1  w3 * w4) - (256 w1  w3) - (128 w2  w3) - (64 w1w1 * w3) - (512 w1  w4) - (256 w2  w4) - (256 w3  w4) - (128 w1w1 * w4) + (64 w1w1 * w3 * w4) - (128 w4w4) - (128 w4w4 * w1) + (32 w4w4 * w1w1) + (192 w4w4  w1 * w2 * w3) + (120 w4w4 * w1w1  w2 * w3) + (64 w4w4 * w1 * w2) + (96 w4w4 * w1w1  w2) - (64 w4w4 * w2) - (64 w4w4 * w3) + (64 w4w4 * w1 * w3) + (96 w4w4 * w1w1  w3) + (32 w4w4 * w2 * w3)
a4 = 128 - (224 w2  w1 * w3) - (56 w2  w1w1  w3) - (448 w2  w1 * w4) - (112 w2  w1w1  w4) - (448 w1  w3 * w4) - (224 w2  w3 * w4) + (640 w1) + (320 w2) + (320 w3) + (640 w4) + (64* w2 * w1) - (112 w2  w1w1) + (32 w1w1) - (224 w2 * w1 * w3 * w4) + (168 w2  w1w1  w3 * w4) + (64 w1  w3) + (32 w2  w3) - (112 w1w1 * w3) + (128 w1  w4) + (64 w2  w4) + (64 w3  w4) - (224 w1w1 * w4) - (112 w1w1 * w3 * w4) + (32 w4w4) - (224 w4w4 * w1) - (56 w4w4 * w1w1) + (168 w4w4  w1 * w2 * w3) + (210 w4w4 * w1w1  w2 * w3) - (112 w4w4 * w1 * w2) + (84 w4w4 * w1w1  w2) - (112 w4w4 * w2) - (112 w4w4 * w3) - (112 w4w4 * w1 * w3) + (84 w4w4 * w1w1  w3) - (56 w4w4 * w2 * w3)
a5 = - (448 w2  w1 * w3 * w4) - (224 w1w1 * w3 * w4) + (384 w3  w4) - (112 w2  w1w1  w3) - (112 w4w4 * w1w1) + (384 w1 * w3) - (224 w4w4 * w1 * w3) + (192 w2  w3) - (224 w2  w1w1  w4) + (192 w1w1) + (252 w4w4 * w1w1  w2 * w3) + (384 w2  w1) - (768) - (224 w4w4 * w1 * w2) - (112 w4w4 * w2 * w3) + (384 w2  w4) + (192 w4w4) + (768 w1  w4)
a6 = 128 + (224 w2  w1 * w3) - (56 w2  w1w1  w3) + (448 w2  w1 * w4) - (112 w2  w1w1  w4) + (448* w1 * w3 * w4) + (224 w2  w3 * w4) - (640 w1) - (320 w2) - (320 w3) - (640 w4) + (64 w2  w1) + (112 w2  w1w1) + (32 w1w1) - (224 w2 * w1 * w3 * w4) - (168 w2  w1w1  w3 * w4) + (64 w1  w3) + (32 w2  w3) + (112 w1w1 * w3) + (128 w1  w4) + (64 w2  w4) + (64 w3  w4) + (224 w1w1 * w4) - (112 w1w1 * w3 * w4) + (32 w4w4) + (224 w4w4 * w1) - (56 w4w4 * w1w1) - (168 w4w4  w1 * w2 * w3) + (210 w4w4 * w1w1  w2 * w3) - (112 w4w4 * w1 * w2) - (84 w4w4 * w1w1  w2) + (112 w4w4 * w2) + (112 w4w4 * w3) - (112 w4w4 * w1 * w3) - (84 w4w4 * w1w1  w3) - (56 w4w4 * w2 * w3)
a7 = 512 + (128 w2  w1 * w3) + (32 w2  w1w1  w3) + (256 w2  w1 * w4) + (64 w2  w1w1  w4) + (256 w1  w3 * w4) + (128 w2  w3 * w4) - (256 w2  w1) + (64 w2  w1w1) - (128 w1w1) + (128 w2 * w1 * w3 * w4) - (192 w2  w1w1  w3 * w4) - (256 w1  w3) - (128 w2  w3) + (64 w1w1 * w3) - (512 w1  w4) - (256 w2  w4) - (256 w3  w4) + (128 w1w1 * w4) + (64 w1w1 * w3 * w4) - (128 w4w4) + (128 w4w4 * w1) + (32 w4w4 * w1w1) - (192 w4w4  w1 * w2 * w3) + (120 w4w4 * w1w1  w2 * w3) + (64 w4w4 * w1 * w2) - (96 w4w4 * w1w1  w2) + (64 w4w4 * w2) + (64 w4w4 * w3) + (64 w4w4 * w1 * w3) - (96 w4w4 * w1w1  w3) + (32 w4w4 * w2 * w3)
a8 = -192 - (48* w2 * w1 * w3) + (52 w2  w1w1  w3) - (96 w2  w1 * w4) + (104 w2  w1w1  w4) - (96 w1  w3 * w4) - (48 w2  w3 * w4) + (320 w1) + (160 w2) + (160* w3) + (320 w4) - (96 w2 * w1) - (24 w2  w1w1) - (48 w1w1) + (208 w2 * w1 * w3 * w4) - (108 w2  w1w1  w3 * w4) - (96 w1  w3) - (48 w2  w3) - (24 w1w1 * w3) - (192* w1 * w4) - (96 w2  w4) - (96* w3 * w4) - (48 w1w1 * w4) + (104 w1w1 * w3 * w4) - (48 w4w4) - (48 w4w4 * w1) + (52 w4w4 * w1w1) - (108 w4w4  w1 * w2 * w3) + (45 w4w4 * w1w1  w2 * w3) + (104 w4w4 * w1 * w2) - (54 w4w4 * w1w1  w2) - (24 w4w4 * w2) - (24 w4w4 * w3) + (104 w4w4 * w1 * w3) - (54 w4w4 * w1w1  w3) + (52 w4w4 * w2 * w3)
a9 = -128 - (64* w2 * w1 * w3) + (24 w2  w1w1  w3) - (128 w2  w1 * w4) + (48 w2  w1w1  w4) - (128 w1  w3 * w4) - (64 w2  w3 * w4) + (64 w2  w1) - (32 w2  w1w1) + (32 w1w1) + (96 w2 * w1 * w3 * w4) - (32 w2  w1w1  w3 * w4) + (64 w1  w3) + (32 w2  w3) - (32 w1w1 * w3) + (128 w1  w4) + (64 w2  w4) + (64 w3  w4) - (64 w1w1 * w4) + (48* w1w1  w3 * w4) + (32 w4w4) - (64 w4w4 * w1) + (24 w4w4 * w1w1) - (32 w4w4  w1 * w2 * w3) + (10 w4w4 * w1w1  w2 * w3) + (48 w4w4 * w1 * w2) - (16 w4w4 * w1w1  w2) - (32 w4w4 * w2) - (32 w4w4 * w3) + (48 w4w4 * w1 * w3) - (16 w4w4 * w1w1  w3) + (24 w4w4 * w2 * w3)
a10 = 64 - (16 w2  w1 * w3) + (4 w2  w1w1  w3) - (32 w2  w1 * w4) + (8 w2  w1w1  w4) - (32 w1  w3 * w4) - (16 w2  w3 * w4) - (64 w1) - (32 w2) - (32 w3) - (64 w4) + (32 w2  w1) - (8 w2  w1w1) + (16 w1w1) + (16 w2 * w1 * w3 * w4) - (4 w2  w1w1  w3 * w4) + (32 w1  w3) + (16 w2  w3) - (8 w1w1 * w3) + (64 w1  w4) + (32 w2  w4) + (32 w3  w4) - (16* w1w1  w4) + (8 w1w1 * w3 * w4) + (16 w4w4) - (16 w4w4 * w1) + (4 w4w4 * w1w1) - (4 w4w4  w1 * w2 * w3) + (w4w4  w1w1  w2 * w3) + (8 w4w4 * w1 * w2) - (2 w4w4 * w1w1  w2) - (8 w4w4 * w2) - (8 w4w4 * w3) + (8 w4w4 * w1 * w3) - (2 w4w4 * w1w1  w3) + (4 w4w4 * w2 * w3)
b0 = 16w4w4
b1 = 32w4w4
b2 = -48w4w4
b3 = -128w4w4
b4 = 32w4w4
b5 = 192w4w4
b6 = 32w4w4
b7 = -128w4w4
b8 = -48w4w4
b9 = 32w4w4
b10 = 16w4w4
Ga = 10^(2/20)
%teste
x=[0 0 0 0 0 0 0 0 0 0 0];
y=[0 0 0 0 0 0 0 0 0 0 0];
t = linspace(0,1,48000);
yy = zeros(1, 48000);
for c= 1:48000
x(1) = sin(2pi100*t(c));
y(1) = (1/a0)(b0x(1) + b1x(2) + b2x(3) + b3x(4) + b4x(5) + b5x(6) + b6x(7) + b7x(8) + b8x(9) + b9x(10) + b10x(11) + a1y(2) + a2y(3) + a3y(4) + a4y(5) + a5y(6) + a6y(7) + a7y(8) + a8y(9) + a9y(10) + a10y(11));
yy(c) = y(1);
  % update x and y data vectors
  for i = 10:-1:1
    x(i+1) = x(i); % store xi
    y(i+1) = y(i); % store yi
  end
end
plot(t,yy)

And nothing but disaster happens when testing with a sine wave with 100Hz:

Signal just gets huge, same happens with other frequencies. What am i doing wrong?

Edit:

sys = tf(10^(2/20).*[b0 b1 b2 b3 b4 b5 b6 b7 b8 b9 b10],[a0 a1 a2 a3 a4 a5 a6 a7 a8 a9 a10],1/fs);
P = pole(sys)
zer=0
zplane(zer,P)
P =
-1.0001 + 0.0001i
  -1.0001 - 0.0001i
  -0.9999 + 0.0001i
  -0.9999 - 0.0001i
   0.9973 + 0.0000i
   0.9973 - 0.0000i
   0.9860 + 0.0000i
   0.9078 + 0.0000i
  -0.0127 + 0.0000i
  -0.0127 - 0.0000i

The poles seem to be almost in the unstable region. Do the poles with real part equal to 1.0001 ruin everything? how can i fix them?

You can't keep dumping code here and expect us to debug it and solve your problem for you. You need to spend time studying filter design, implementation and properties with maybe simpler examples and debugging your own code thoroughly. Use @Hilmar's recommendation on your previous question as a general guideline. — Jdip, Dec 12 '22 at 16:29
If you're not asked to use that high order system showed in your post, there already are few digital A-weighting filter implementation threads on this forum. Example. My implementation there gives good match against the definition of A-weighting filter but, you can also edit the code to get exactly same response seen in your linked plot (just use Butterworth/BLT instead of MIM for the LPF). — Juha P, Dec 13 '22 at 08:50
We've been over this in your other question https://dsp.stackexchange.com/questions/85729/trying-to-implement-a-digital-lpf. Your filter conversion is clearly wrong. The poles are outside the unit circle and hence the filter is unstable. You need to debug this one step at a time as explained in the answer to your other question. — Hilmar, Dec 13 '22 at 09:05

TimWescott · Accepted Answer · 2022-12-15T00:08:49.157

IIR filters should, almost exactly universally* be broken down into first- and second-order sections and cascaded**. The sensitivity of filter pole locations to coefficient values goes up with filter order; even the slightest rounding error will screw up a 10th-order filter.
If you have the signal processing toolbox, you should use the built-in IIR filter function.
1. If you don't, you should still vectorize a and b

My second-choice recommendation is to vectorize a and b and use a polynomial root-finder to verify that the poles and zeros are in sensible locations, then try to find the bug that's making them wrong.

The poles seem to be almost in the unstable region.

No, the poles are in the unstable region. Anything outside the unit circle ($|z| > 1$) is unstable.

Do the poles with real part equal to 1.0001 ruin everything?

No. They are diagnostic of the fact that everything is ruined. You have some underlying trouble that needs to be fixed.

how can i fix them?

Sensibly, take my first-choice recommendation, below, or do what I'd be tempted to do myself. That's why I'm recommending them.
- Note that someone has already volunteered you a link to a solution.
If you must persist, try a math package that has arbitrarily high precision (there may be a Matlab extension that does this), and try the root-finding at higher precision. But be aware that you're going down a rabbit-hole, and in a world where people publish designs for this sort of thing for free and because it's fun, it's a pointless rabbit-hole.

My first-choice recommendation is to try to find a paper that gives you the poles and zeros of a cascade of 2nd-order filters.

What I'd be tempted to do if my first choice didn't work out would be to fit my own IIR filters to the recommended filter shape. Even if I did have that first-choice reference, I may do it anyway and compare which looks better.

* unless you're an absolute freaking expert and willing to argue with your colleagues and expect to win, you should read this as "absolutely universally".

** or, rarely, cascade-parallel -- this would apply if you have a wide notch filter or similar.

Agree. The MATLAB, Octave and Python Scipy.signal tools support "sos" mode for IIR filter constructions: returning the coeffiicents as second order sections with a reasonable grouping. I believe it comes down to selecting the complex conjugate pair of poles that is closest in Euclidean distance on the z-plane to the complex conjugate pair of zeros (minimizing the gain for that section and ultimately dynamic range requirements for the filter implementation, when combined with proper scaling between sections). — Dan Boschen, Dec 12 '22 at 17:54
Can you see my new post edit? Perhaps the poles with real part equal to 1.0001 ruin everything, how can i fix them? — Scipio, Dec 14 '22 at 23:36

Trying to implement a digital A frequency filter

1 Answers1

Linked