Trying to implement a digital A frequency filter (continuation)

Question

Following this post some user advised that the filter developed there should be broken down into first and second-order sections and cascaded. In a paper i found the same thing is told:

So i made what the authors told in matlab:

B = [1 0 ; 1 0.3 ; 1 -1; 1 -1 ; 1 -1; 1 -1];
A = [1 -0.2025 ; 1 -0.2025 ; 1 -0.9860 ; 1 -0.9079 ; 1 -0.9973 ; 1 -0.9973];

[H1,f1] = freqz([B(1,1) B(1,2)] ,[A(1,1) A(1,2)],512,48000);
[H2,f2] = freqz([B(2,1) B(2,2)] ,[A(2,1) A(2,2)],512,48000);
[H3,f3] = freqz([B(3,1) B(3,2)] ,[A(3,1) A(3,2)],512,48000);
[H4,f4] = freqz([B(4,1) B(4,2)] ,[A(4,1) A(4,2)],512,48000);
[H5,f5] = freqz([B(5,1) B(5,2)] ,[A(5,1) A(5,2)],512,48000);
[H6,f6] = freqz([B(6,1) B(6,2)] ,[A(6,1) A(6,2)],512,48000);

mag = abs( H1.*H2.*H3.*H4.*H5.*H6);

magdb = 20*log(mag);
semilogx(f2,magdb)

ylim([-70 40])
title('Digital implementation of the A-weighting filter (fs = 48 kHz)')
xlabel('Frequency (Hz)')
ylabel('Gain(dB)')

And obtained this:

Something is wrong at the centre frequencies because the gain is 10dB when it should be close to 0dB like this:

Why is that happening?

In the other post 2 poles where outside the unit circle, here all poles are inside.

Code developed:

sys1 = tf([B(1,1) B(1,2)] ,[A(1,1) A(1,2)]);
sys2 = tf([B(2,1) B(2,2)] ,[A(2,1) A(2,2)]);
sys3 = tf([B(3,1) B(3,2)] ,[A(3,1) A(3,2)]);
sys4 = tf([B(4,1) B(4,2)] ,[A(4,1) A(4,2)]);
sys5 = tf([B(5,1) B(5,2)] ,[A(5,1) A(5,2)]);
sys6 = tf([B(6,1) B(6,2)] ,[A(6,1) A(6,2)]);



P = pole (sys1)
P = pole (sys2)
P = pole (sys3)
P = pole (sys4)
P = pole (sys5)
P = pole (sys6)

Output:

P = pole (sys1)
P = pole (sys2)
P = pole (sys3)
P = pole (sys4)
P = pole (sys5)
P = pole (sys6)

P = 0.2025
P = 0.2025
P = 0.9860
P = 0.9079
P = 0.9973
P = 0.9973

So i tried to test just the first order section to see where it goes:

figure()
mag = abs( H1);
magdb = 20*log(mag);
semilogx(f1,magdb);

ylim([-70 40])
xlabel('Frequency (Hz)')
ylabel('Gain(dB)')
grid on

I tested it with a sine wave with amplitude one where i will change its frequency:

x=[0 0];
y=[0 0];

t = linspace(0,1,fs);

yy = zeros(1,fs);


for c= 1:fs
x(1) = sin(2*pi*10000*t(c));


y(1) = (1/A(1,1))*(B(1,1)*x(1) + B(1,2)*x(2) + A(1,2)*y(2) ); 

yy(c) = y(1);
  % update x and y data vectors
  for i = 2:-1:1
    x(i+1) = x(i); % store xi
    y(i+1) = y(i); % store yi
  end
  
 
end

figure()
plot(t,yy)

For 100Hz i got this:

The gain for that frequency is around 5dB so i dont know why the amplitude is less than 1.

For 1000Hz i got this:

Same as before.

For 6500Hz:

The amplotitude should drop by now but it seems to get a bit higher

For 10000Hz:

For some reason the amplitude is getting higher.

15000Hz:

It gets higher, this section doesn't seem to behave like its transfer function... Why?

Also, now i have 6 cascaded filters, what's the procedure to make this for all 6 of them at once?

Answer 1

I don't know that this is the issue, but the original transfer function shown has a 0.9097 in the denominator which was transcribed incorrectly in the first order sections to 0.9079. (One or the other is in error). I would also be very suspicious of round-off error when only 4 significant places are shown for the filter coefficients.

I confirmed different results for 100 Hz so perhaps an error in the generation of the test waveform or the implementation of the filter. Rather than debugging that, I suggest using the filter command available in MATLAB as well as MATLAB's vector processing.

Here is the test for 100 Hz along with the code using the approach I would take:

fs = 48000;
nsamps = 48000;
t = [0: nsamps-1]/fs;
fout = 100;
x = sin(2*pi*fout*t);

b = [1 0];
a = [1 -0.2025];
y = filter(b, a, x);

figure
plot(y)

With the following result: