Why does the filter function in matlab get different results from the IIR filters using direct form?

Question

I have 6 first order filters, each one has 2 coefficients (B(x1,x1), B(x2,x2) for the numerator and A(y1,y1),A(y2,y2) for the denominator).

B = [0.6 0 ; 1 0.3 ; 1 -1; 1 -1 ; 1 -1; 1 -1]; % b(1,1) ajustado para 0.6
A = [1 -0.2025 ; 1 -0.2025 ; 1 -0.9860 ; 1 -0.9097 ; 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); % pus um ajuste de 0.5

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

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

This is the transfer function i obtain (ignore the title):

And here are all the poles (each pole corresponds to one first order filter):

P =

   0.202500000000000


P =

   0.202500000000000


P =

   0.986000000000000


P =

   0.909700000000000


P =

   0.997300000000000


P =

   0.997300000000000

Every pole is inside the unit circle.

This is the code and frequency response of the third first order filter:

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

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

Using the filter function i obtain the following for a 10000Hz sin wave:

t = linspace(0,1,fs);
x = sin(2*pi*10000*t);

y = filter([B(3,1) B(3,2)],[A(3,1) A(3,2)],x);

figure()
plot(t,x)
hold on
plot(t,y)
legend('Input Data','Filtered Data')

This is exactly what i'm expecting.

But now if i use the IIR Filter coeficcent using the past values i get something different:

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(3,1))*(B(3,1)*x(1) + B(3,2)*x(2) + A(3,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)

Why does this happen? How can i get the same response as the filter function?

Answer 1

I think the y(1) equation should be

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

then I get

Note also, the update should just be

x(2) = x(1)
y(2) = y(1)

otherwise the x and y matrices grow to be $3 \times 1$ after the first time through the inner for loop.

Answer 2

Both your formula and the other answer formula for y(1) is wrong. Here is the correct formula:

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

Relative to the formula in your question, notice the subtraction rather than addition of the final term.

An alternative way to implement it that matches Matlab's filter function can be found in the algorithm section of the corresponding documentation.