C equivalent to MATLAB “filter”?

Why is my MATLAB "filter" giving different result than equivalent C code?

MATLAB:

[B_coeffs, A_coeffs ] = butter(4, 100/(16000/2), 'high');

state  = zeros( 4, 1 );
input = zeros( 64,1 );

for i=1:64
    input(i)=i;
end

[filtered_output, state] = filter( B_coeffs, A_coeffs, input, state );

C:

int main(...)
{
    for(int test=0; test<64;test++)
        Xin[test]=test+1;
    ...
    high_pass_filter_init(...)
    high_pass_filter_do(...)
}

// Do the filtering
void high_pass_filter_do( t_high_pass_filter* hpf, float *Xin, float *Yout )
{
    double Xi, Yi;

      double z0 = hpf->state[0],
             z1 = hpf->state[1],
             z2 = hpf->state[2],
             z3 = hpf->state[3];

      unsigned int N = 64;
      unsigned int i = 0;

      for(i=0;i<N;i++)
      {
          Xi = Xin[i];

          Yi = hpf->B_coeffs[0] * Xi + z0;
          z0 = hpf->B_coeffs[1] * Xi + z1 - hpf->A_coeffs[1] * Yi;
          z1 = hpf->B_coeffs[2] * Xi + z2 - hpf->A_coeffs[2] * Yi;
          z2 = hpf->B_coeffs[3] * Xi + z3 - hpf->A_coeffs[3] * Yi;
          z3 = hpf->B_coeffs[4] * Xi      - hpf->A_coeffs[4] * Yi;

          Yout[i] = (float) Yi;
      }

      hpf->state[0] = z0;
      hpf->state[1] = z1;
      hpf->state[2] = z2;
      hpf->state[3] = z3;

      return;
  }

// Filter struct
typedef struct
{
    float A_coeffs[5];
    float B_coeffs[5];
    float state[4];    
} t_high_pass_filter;

// Create filter
void high_pass_filter_init( t_high_pass_filter* hpf)
{

      hpf->A_coeffs[0] =  1.0000;
      hpf->A_coeffs[1] = -3.8974;
      hpf->A_coeffs[2] =  5.6974;
      hpf->A_coeffs[3] = -3.7025;
      hpf->A_coeffs[4] =  0.9025;

      hpf->B_coeffs[0] =  0.9500;
      hpf->B_coeffs[1] = -3.7999;
      hpf->B_coeffs[2] =  5.6999;
      hpf->B_coeffs[3] = -3.7999;
      hpf->B_coeffs[4] =  0.9500;

      hpf->state[0] = 0.0;
      hpf->state[1] = 0.0;
      hpf->state[2] = 0.0;
      hpf->state[3] = 0.0;
  }

The outputs are similar at the beggining, but then they diverge:

MATLAB:        C: 
----------------------------
0.9500         0.9500
1.8025         1.8026
2.5625         2.5631
3.2350         3.2369
3.8247         3.8292
4.3360         4.3460
4.7736         4.7930
5.1416         5.1767
5.4442         5.5035
5.6854         5.7807
5.8691         6.0156
5.9991         6.2162
6.0788         6.3909
6.1119         6.5487
6.1016         6.6989
6.0511         6.8514
5.9637         7.0167
5.8421         7.2057
5.6894         7.4298
5.5083         7.7009
5.3013         8.0314
5.0710         8.4342
4.8199         8.9225
4.5501         9.5101
4.2640        10.2110
3.9637        11.0399
3.6511        12.0115
3.3281        13.1412
2.9965        14.4443
2.6582        15.9368
2.3146        17.6347
1.9674        19.5543
1.6180        21.7122
1.2677        24.1250
0.9179        26.8095
0.5698        29.7829
0.2245        33.0621
-0.1169        36.6643
-0.4535        40.6067
-0.7842        44.9066
-1.1084        49.5812
-1.4251        54.6477
-1.7336        60.1232
-2.0333        66.0249
-2.3236        72.3697
-2.6039        79.1746
-2.8738        86.4562
-3.1327        94.2313
-3.3804       102.5161
-3.6164       111.3270
-3.8405       120.6801
-4.0524       130.5911
-4.2520       141.0756
-4.4390       152.1490
-4.6134       163.8264
-4.7750       176.1225
-4.9239       189.0520
-5.0600       202.6291
-5.1832       216.8677
-5.2938       231.7814
-5.3917       247.3837
-5.4771       263.6875
-5.5501       280.7055
-5.6108       298.4500

What am I doing wrong??

Answer 1

[Jan]

I get nearly the same output in Matlab, if I use the same coefficients as declared in your C-code:

A = [1.0000; -3.8974; 5.6974; -3.7025; 0.9025];
B = [0.9500; -3.7999; 5.6999; -3.7999; 0.9500];
filter(B, A, 1:64)

Note that the parameters replied by butter are (about - I do not have the SPT installed and use my own implementation of Butter with tiny differences in the output):

A = [1.0, -3.89738598246839, ...
     5.69739000394907, -3.70246714846478, 0.902465387346562]
B = [0.9499817826393, -3.7999271305572, ...
     5.6998906958358, -3.7999271305572, 0.9499817826393]

And due to the different parameters, you get a different result.

[Danijel]

It's so sensitive to so small coefficient differences to produce totaly different result??

[Jan]

The differences are in the 5th digit after the point. Because you sum over 5 elements, loosing a valid digit after each 4 steps is plausible, becasue the differences accumulate. This sensitivity is in the expected order.