Trying to implement a digital A frequency filter

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 = 2*tan(pi*(f1/fs));
w2 = 2*tan(pi*(f2/fs));
w3 = 2*tan(pi*(f3/fs));
w4 = 2*tan(pi*(f4/fs));

%testeb2 = 2*tan(pi*(250/1000))*2*tan(pi*(250/1000))*(1/sqrt(2))
% Filter coefficients for the A weighting filter

a0 = 64 + (16*w2*w1*w3) + (4*w2*w1*w1*w3) + (32*w2*w1*w4) + (8*w2*w1*w1*w4) + (32* w1 * w3 * w4) + (16 *w2 * w3 * w4) + (64 *w1) + (32 *w2) + (32 *w3) + (64 *w4) + (32 *w2 * w1) + (8 *w2 * w1*w1) + (16 *w1*w1) + (16 *w2 * w1 * w3 * w4) + (4 *w2 * w1*w1 * w3 * w4) + (32 *w1 * w3) + (16 *w2 * w3) + (8 *w1*w1 * w3) + (64 *w1 * w4) + (32 *w2 * w4) + (32 *w3 * w4) + (16 *w1*w1 * w4) + (8 *w1*w1 * w3 * w4) + (16 *w4*w4) + (16 *w4*w4 * w1) + (4 *w4*w4 * w1*w1) + (4 *w4*w4 * w1 * w2 * w3) + (w4*w4 * w1*w1 * w2 * w3) + (8 *w4*w4 * w1 * w2) + (2 *w4*w4 * w1*w1 * w2) + (8 *w4*w4 * w2) + (8 *w4*w4 * w3) + (8 *w4*w4 * w1 * w3) + (2 *w4*w4 * w1*w1 * w3) + (4 *w4*w4 * w2 * w3)

a1 = -128 + (64 *w2 * w1 * w3) + (24* w2 * w1*w1 * w3) + (128 *w2 * w1 * w4) + (48 *w2 * w1*w1 * w4) + (128 *w1 * w3 * w4) + (64 *w2 * w3 * w4) + (64 *w2 * w1) + (32 *w2 * w1*w1) + (32 *w1*w1) + (96 *w2 * w1 * w3 * w4) + (32 *w2 * w1*w1 * w3 * w4) + (64 *w1 * w3) + (32 *w2 * w3) + (32 *w1*w1 * w3) + (128 *w1 * w4) + (64 *w2 * w4) + (64 *w3 * w4) + (64 *w1*w1 * w4) + (48 *w1*w1 * w3 * w4) + (32 *w4*w4) + (64 *w4*w4 * w1) + (24 *w4*w4 * w1*w1) + (32 *w4*w4 * w1 * w2 * w3) + (10 *w4*w4 * w1*w1 * w2 * w3) + (48 *w4*w4 * w1 * w2) + (16 *w4*w4 * w1*w1 * w2) + (32 *w4*w4 * w2) + (32 *w4*w4 * w3) + (48 *w4*w4 * w1 * w3) + (16 *w4*w4 * w1*w1 * w3) + (24 *w4*w4 * w2 * w3)

a2 = -192 + (48 *w2 * w1 * w3) + (52 *w2 * w1*w1 * w3) + (96 *w2 * w1 * w4) + (104 *w2 * w1*w1 * w4) + (96 *w1 * w3 * w4) + (48 *w2 * w3 * w4) - (320 *w1) - (160 *w2) - (160 *w3) - (320 *w4) - (96 *w2 * w1) + (24 *w2 * w1*w1) - (48 *w1*w1) + (208 *w2 * w1 * w3 * w4) + (108 *w2 * w1*w1 * w3 * w4) - (96 *w1 * w3) - (48 *w2 * w3) + (24 *w1*w1 * w3) - (192 *w1 * w4) - (96 *w2 * w4) - (96 *w3 * w4) + (48 *w1*w1 * w4) + (104 *w1*w1 * w3 * w4) - (48 *w4*w4) + (48 *w4*w4 * w1) + (52 *w4*w4 * w1*w1) + (108 *w4*w4 * w1 * w2 * w3) + (45 *w4*w4 * w1*w1 * w2 * w3) + (104 *w4*w4 * w1 * w2) + (54 *w4*w4 * w1*w1 * w2) + (24 *w4*w4 * w2) + (24 *w4*w4 * w3) + (104 *w4*w4 * w1 * w3) + (54 *w4*w4 * w1*w1 * w3) + (52 *w4*w4 * w2 * w3)

a3 = 512 - (128 *w2 * w1 * w3) + (32 *w2 * w1*w1 * w3) - (256 *w2 * w1 * w4) + (64 *w2 * w1*w1 * w4) - (256 *w1 * w3 * w4) - (128 *w2 * w3 * w4) - (256 *w2 * w1) - (64 *w2 * w1*w1) - (128 *w1*w1) + (128 *w2 * w1 * w3 * w4) + (192* w2 * w1*w1 * w3 * w4) - (256 *w1 * w3) - (128 *w2 * w3) - (64 *w1*w1 * w3) - (512 *w1 * w4) - (256 *w2 * w4) - (256 *w3 * w4) - (128 *w1*w1 * w4) + (64 *w1*w1 * w3 * w4) - (128 *w4*w4) - (128 *w4*w4 * w1) + (32 *w4*w4 * w1*w1) + (192 *w4*w4 * w1 * w2 * w3) + (120 *w4*w4 * w1*w1 * w2 * w3) + (64 *w4*w4 * w1 * w2) + (96 *w4*w4 * w1*w1 * w2) - (64 *w4*w4 * w2) - (64 *w4*w4 * w3) + (64 *w4*w4 * w1 * w3) + (96 *w4*w4 * w1*w1 * w3) + (32 *w4*w4 * w2 * w3)

a4 = 128 - (224 *w2 * w1 * w3) - (56 *w2 * w1*w1 * w3) - (448 *w2 * w1 * w4) - (112 *w2 * w1*w1 * w4) - (448 *w1 * w3 * w4) - (224 *w2 * w3 * w4) + (640 *w1) + (320 *w2) + (320 *w3) + (640 *w4) + (64* w2 * w1) - (112 *w2 * w1*w1) + (32 *w1*w1) - (224 *w2 * w1 * w3 * w4) + (168 *w2 * w1*w1 * w3 * w4) + (64 *w1 * w3) + (32 *w2 * w3) - (112 *w1*w1 * w3) + (128 *w1 * w4) + (64 *w2 * w4) + (64 *w3 * w4) - (224 *w1*w1 * w4) - (112 *w1*w1 * w3 * w4) + (32 *w4*w4) - (224 *w4*w4 * w1) - (56 *w4*w4 * w1*w1) + (168 *w4*w4 * w1 * w2 * w3) + (210 *w4*w4 * w1*w1 * w2 * w3) - (112 *w4*w4 * w1 * w2) + (84 *w4*w4 * w1*w1 * w2) - (112 *w4*w4 * w2) - (112 *w4*w4 * w3) - (112 *w4*w4 * w1 * w3) + (84 *w4*w4 * w1*w1 * w3) - (56 *w4*w4 * w2 * w3)

a5 = - (448 *w2 * w1 * w3 * w4) - (224 *w1*w1 * w3 * w4) + (384 *w3 * w4) - (112 *w2 * w1*w1 * w3) - (112 *w4*w4 * w1*w1) + (384 *w1 * w3) - (224 *w4*w4 * w1 * w3) + (192 *w2 * w3) - (224 *w2 * w1*w1 * w4) + (192 *w1*w1) + (252 *w4*w4 * w1*w1 * w2 * w3) + (384 *w2 * w1) - (768) - (224 *w4*w4 * w1 * w2) - (112 *w4*w4 * w2 * w3) + (384 *w2 * w4) + (192 *w4*w4) + (768 *w1 * w4)

a6 = 128 + (224 *w2 * w1 * w3) - (56 *w2 * w1*w1 * w3) + (448 *w2 * w1 * w4) - (112 *w2 * w1*w1 * w4) + (448* w1 * w3 * w4) + (224 *w2 * w3 * w4) - (640 *w1) - (320 *w2) - (320 *w3) - (640 *w4) + (64 *w2 * w1) + (112 *w2 * w1*w1) + (32 *w1*w1) - (224 *w2 * w1 * w3 * w4) - (168 *w2 * w1*w1 * w3 * w4) + (64 *w1 * w3) + (32 *w2 * w3) + (112 *w1*w1 * w3) + (128 *w1 * w4) + (64 *w2 * w4) + (64 *w3 * w4) + (224 *w1*w1 * w4) - (112 *w1*w1 * w3 * w4) + (32 *w4*w4) + (224 *w4*w4 * w1) - (56 *w4*w4 * w1*w1) - (168 *w4*w4 * w1 * w2 * w3) + (210 *w4*w4 * w1*w1 * w2 * w3) - (112 *w4*w4 * w1 * w2) - (84 *w4*w4 * w1*w1 * w2) + (112 *w4*w4 * w2) + (112 *w4*w4 * w3) - (112 *w4*w4 * w1 * w3) - (84 *w4*w4 * w1*w1 * w3) - (56 *w4*w4 * w2 * w3)

a7 = 512 + (128 *w2 * w1 * w3) + (32 *w2 * w1*w1 * w3) + (256 *w2 * w1 * w4) + (64 *w2 * w1*w1 * w4) + (256 *w1 * w3 * w4) + (128 *w2 * w3 * w4) - (256 *w2 * w1) + (64 *w2 * w1*w1) - (128 *w1*w1) + (128 *w2 * w1 * w3 * w4) - (192 *w2 * w1*w1 * w3 * w4) - (256 *w1 * w3) - (128 *w2 * w3) + (64 *w1*w1 * w3) - (512 *w1 * w4) - (256 *w2 * w4) - (256 *w3 * w4) + (128 *w1*w1 * w4) + (64 *w1*w1 * w3 * w4) - (128 *w4*w4) + (128 *w4*w4 * w1) + (32 *w4*w4 * w1*w1) - (192 *w4*w4 * w1 * w2 * w3) + (120 *w4*w4 * w1*w1 * w2 * w3) + (64 *w4*w4 * w1 * w2) - (96 *w4*w4 * w1*w1 * w2) + (64 *w4*w4 * w2) + (64 *w4*w4 * w3) + (64 *w4*w4 * w1 * w3) - (96 *w4*w4 * w1*w1 * w3) + (32 *w4*w4 * w2 * w3)

a8 = -192 - (48* w2 * w1 * w3) + (52 *w2 * w1*w1 * w3) - (96 *w2 * w1 * w4) + (104 *w2 * w1*w1 * w4) - (96 *w1 * w3 * w4) - (48 *w2 * w3 * w4) + (320 *w1) + (160 *w2) + (160* w3) + (320 *w4) - (96 *w2 * w1) - (24 *w2 * w1*w1) - (48 *w1*w1) + (208 *w2 * w1 * w3 * w4) - (108 *w2 * w1*w1 * w3 * w4) - (96 *w1 * w3) - (48 *w2 * w3) - (24 *w1*w1 * w3) - (192* w1 * w4) - (96 *w2 * w4) - (96* w3 * w4) - (48 *w1*w1 * w4) + (104 *w1*w1 * w3 * w4) - (48 *w4*w4) - (48 *w4*w4 * w1) + (52 *w4*w4 * w1*w1) - (108 *w4*w4 * w1 * w2 * w3) + (45 *w4*w4 * w1*w1 * w2 * w3) + (104 *w4*w4 * w1 * w2) - (54 *w4*w4 * w1*w1 * w2) - (24 *w4*w4 * w2) - (24 *w4*w4 * w3) + (104 *w4*w4 * w1 * w3) - (54 *w4*w4 * w1*w1 * w3) + (52 *w4*w4 * w2 * w3)

a9 = -128 - (64* w2 * w1 * w3) + (24 *w2 * w1*w1 * w3) - (128 *w2 * w1 * w4) + (48 *w2 * w1*w1 * w4) - (128 *w1 * w3 * w4) - (64 *w2 * w3 * w4) + (64 *w2 * w1) - (32 *w2 * w1*w1) + (32 *w1*w1) + (96 *w2 * w1 * w3 * w4) - (32 *w2 * w1*w1 * w3 * w4) + (64 *w1 * w3) + (32 *w2 * w3) - (32 *w1*w1 * w3) + (128 *w1 * w4) + (64 *w2 * w4) + (64 *w3 * w4) - (64 *w1*w1 * w4) + (48* w1*w1 * w3 * w4) + (32 *w4*w4) - (64 *w4*w4 * w1) + (24 *w4*w4 * w1*w1) - (32 *w4*w4 * w1 * w2 * w3) + (10 *w4*w4 * w1*w1 * w2 * w3) + (48 *w4*w4 * w1 * w2) - (16 *w4*w4 * w1*w1 * w2) - (32 *w4*w4 * w2) - (32 *w4*w4 * w3) + (48 *w4*w4 * w1 * w3) - (16 *w4*w4 * w1*w1 * w3) + (24 *w4*w4 * w2 * w3)

a10 = 64 - (16 *w2 * w1 * w3) + (4 *w2 * w1*w1 * w3) - (32 *w2 * w1 * w4) + (8 *w2 * w1*w1 * w4) - (32 *w1 * w3 * w4) - (16 *w2 * w3 * w4) - (64 *w1) - (32 *w2) - (32 *w3) - (64 *w4) + (32 *w2 * w1) - (8 *w2 * w1*w1) + (16 *w1*w1) + (16 *w2 * w1 * w3 * w4) - (4 *w2 * w1*w1 * w3 * w4) + (32 *w1 * w3) + (16 *w2 * w3) - (8 *w1*w1 * w3) + (64 *w1 * w4) + (32 *w2 * w4) + (32 *w3 * w4) - (16* w1*w1 * w4) + (8 *w1*w1 * w3 * w4) + (16 *w4*w4) - (16 *w4*w4 * w1) + (4 *w4*w4 * w1*w1) - (4 *w4*w4 * w1 * w2 * w3) + (w4*w4 * w1*w1 * w2 * w3) + (8 *w4*w4 * w1 * w2) - (2 *w4*w4 * w1*w1 * w2) - (8 *w4*w4 * w2) - (8 *w4*w4 * w3) + (8 *w4*w4 * w1 * w3) - (2 *w4*w4 * w1*w1 * w3) + (4 *w4*w4 * w2 * w3)

b0 = 16*w4*w4

b1 = 32*w4*w4

b2 = -48*w4*w4

b3 = -128*w4*w4

b4 = 32*w4*w4

b5 = 192*w4*w4

b6 = 32*w4*w4

b7 = -128*w4*w4

b8 = -48*w4*w4

b9 = 32*w4*w4

b10 = 16*w4*w4

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(2*pi*100*t(c));

y(1) = (1/a0)*(b0*x(1) + b1*x(2) + b2*x(3) + b3*x(4) + b4*x(5) + b5*x(6) + b6*x(7) + b7*x(8) + b8*x(9) + b9*x(10) + b10*x(11) + a1*y(2) + a2*y(3) + a3*y(4) + a4*y(5) + a5*y(6) + a6*y(7) + a7*y(8) + a8*y(9) + a9*y(10) + a10*y(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
Commented Dec 12, 2022 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). Commented Dec 13, 2022 at 8:50
• We've been over this in your other question dsp.stackexchange.com/questions/85729/…. 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. Commented Dec 13, 2022 at 9:05

1. 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.
2. 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). Commented Dec 12, 2022 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? Commented Dec 14, 2022 at 23:36