0
$\begingroup$

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

enter image description here

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

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

enter image description here enter image description here

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:

enter image description here

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

enter image description here

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?

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ 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. $\endgroup$
    – Jdip
    Commented Dec 12, 2022 at 16:29
  • $\begingroup$ 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). $\endgroup$
    – Juha P
    Commented Dec 13, 2022 at 8:50
  • $\begingroup$ 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. $\endgroup$
    – Hilmar
    Commented Dec 13, 2022 at 9:05

1 Answer 1

Reset to default
3
$\begingroup$
  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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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). $\endgroup$
    – Dan Boschen
    Commented Dec 12, 2022 at 17:54
  • $\begingroup$ Can you see my new post edit? Perhaps the poles with real part equal to 1.0001 ruin everything, how can i fix them? $\endgroup$
    – Scipio
    Commented Dec 14, 2022 at 23:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.